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The  ground-state  structural  energies  of  several  binary  compounds 
in  the  systems  zirconium-niobium,  niobium-ruthenium,  and 
zirconium-ruthenium  were  calculated  using  first- principles  total- 
energy  electronic-structure  calculations  with  the  linearized  muffin-tin 
orbital  (LMTO)  approximation.  These  energies  were  used  in  a  cluster 
expansion  to  determine  the  ground-state  configurations  of  each  of  the 
three  binary  systems.  Notably,  the  calculations  imply  that  the  stable 
equilibrium  structure  of  the  compound  NbRu3  is  DO  19.  The 
calculated  excess  structural  energies  are  in  good  agreement  with 
reported  measurements.  The  long-term  objective  of  this  approach  is 
the  calculation  of  equilibrium  thermomechanical  properties  of  alloys 
without  resort  to  experimental  measurements,  for  example,  phase 
diagrams.  Finite  temperature  phase  stability  calculations  are 
implemented  using  a  generalized  I  sing  model  and  the  cluster  variation 
method  (CVM).  Vibrational  free  energy  is  included  using  the  Debye- 
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Gruneisen  approximation.  Predicted  phase  equilibria  are  in  excellent 
agreement  with  measurements. 
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INTRODUCTION 


Progress  in  a  long-standing  objective  in  computational 
materials  science,  the  calculation  of  alloy  phase  diagrams  from  first- 
principles.  has  been  paced  by  the  associated  problems  of  accurate  yet 
efficient  total-energy  calculations  and  statistical  approximations  of 
the  alloy  partition  function.  The  former  is  well  understood,  and  the 
evolution  of  Density-Functional  Theory  (DFT)  [11  to  computationally 
suitable  methods  using  the  Local-Density  Approximation  (LDA)  [2]  is 
a  renowned  (though,  arguably,  a  qualified)  success.  The  high 
precision  required  of  the  total-energy  band-structure  methods  for 
phase-equilibrium  calculations  is  attainable  for  most  chemical 
systems.  The  connection  between  the  total-energy  calculations  and  a 
general  solution  to  the  Ising  problem  is  provided  by  a  cluster 
expansion  in  a  series  of  multi-atom  interaction  energies  [3].  This 
rigorous  expansion  may  be  inverted  to  calculate  the  total  energy  of 
any  lattice  configuration  using  the  multi-atom  interactions.  Two 
main  aspects  of  the  cluster  expansion  method  remain  unresolved:  the 
convergence  of  the  expansion  and  the  effect  of  strain  energy  on  the 
convergence.  These  two  issues  are  central  to  this  thesfr .  and  a  major 
portion  of  this  study  is  devoted  toward  a  resolution. 

First-principles  phase  diagrams  have  been  computed  with 
partial  success  for  noble  metals  (3-5),  semiconductors  [6],  metal-oxide 
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superconductors  (7],  ordering  systems  [8-17],  and  systems  dominated 
by  the  effects  of  coherency  [18]  or  surfaces  [19].  Yet  important 
aspects  are  missing  from  each  of  these  studies.  For  example, 
vibrational  free  energy  is  rarely  included,  except  using  strictly 
empirical  means.  Incorporation  of  elastic  relaxation,  critical  in  some 
cases,  is  far  from  satisfactory  resolution.  Local-volume  relaxation 
has  been  included  by  somewhat  different  methods,  as  have 
distortional  modes,  but  as  yet  even  this  seemingly  simple  problem  has 
no  rigorous  treatment.  For  even  the  simplest  systems,  successful 
phase-diagram  calculation  requires  much  experience  to  complement 
the  impressive  set  of  computational  tools  now  available. 

In  this  study  phase  stability  in  the  systems  zirconium- 
niobium,  niobium-ruthenium,  and  zirconium-ruthenium  are 
investigated  via  the  cluster  expansion  using  total  energies  obtained 
with  the  Linear-Muffln-Tin-Orbital  method  [20]  in  the  Atomic- 
Sphere  Approximation  (LMTO-ASA).  The  systems  are  well  chosen  to 
test  at  different  levels  the  method  for  calculating  phase  equilibrium 
and  for  comparison  with  experimental  data  (though  few  data  have 
been  reported  for  Nb-Ru).  The  Zr-Ru  system  is  simple  with  only  one 
intermediate  ground-state  ordered  compound  —  ZrRu  with  the  B2 
structure.  Zr-Nb  is  also  a  simple  system,  with  a  broad  miscibility  gap, 
and  is  a  good  test  for  the  accuracy  of  the  method.  The  phase 
transition  of  Zr  from  hep  to  bcc  requires  some  treatment  of 
vibrational  free  energy.  The  Nb-Ru  system  is  complex,  and  several 
subtle  aspects  of  phase  equilibrium  have  eluded  experimental 
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Investigators.  The  cluster  expansion  of  the  formation  energy  of  alloys 
In  this  system  are  used  to  resolve  some  of  these  controversial  aspects. 

The  Ising  Hamiltonian  for  a  given  lattice  configuration  can 
be  expressed  through  a  series  expansion  in  the  cluster  correlation 
functions  (variables  describing  the  geometry  in  an  average  sense)  of 
each  type  of  cluster  of  lattice  sites  comprising  the  lattice  (crystal). 
This  expansion  is  rigorous  and  exact  in  the  thermodynamical  limit 
(21,  22).  The  coefficients  of  the  cluster  expansion,  known  as  the 
effective  chemical  interactions  provide  the  connection  between  the 
first-principles  total-energy  calculations  and  the  approximate 
expansion  solutions  to  the  Ising  model.  These  interaction  energies 
can  be  obtained  through  inversion  of  the  cluster  expansions  of  the 
total  energies  of  several  perfectly  ordered  structures.  The  application 
of  the  cluster  expansion  using  ab  initio  total-energy  calculations  was 
first  made  by  Connolly  and  Williams  [23],  and  is  now  a  standard 
method  in  first-principles  phase-equilibrium  calculations.  For  a 
perfectly  ordered  structure  the  correlations  may  be  ascertained  a 
priori  Thus,  using  the  interaction  energies,  the  cluster  expansion 
may  be  used  to  ascertain  the  formation  energy  of  any  ordered 
structure.  This  obviates  total-energy  calculations  of  every  likely 
lattice  configuration  in  the  process  of  determining  the  stable- 
equilibrium  ground-state  structures. 

The  cluster-expansion  method  as  used  by  Connolly  and 
Williams  includes  no  implicit  prescription  for  the  volume  dependence 
of  the  formation  energy  of  the  partially  ordered  alloy.  Sever 
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techniques  for  including  elastic  relaxation  into  the  minimization  of 
the  configurational  energy  have  been  studied.  The  so-called  "G-e" 
method  [6,  10]  is  based  upon  the  stringent  presupposition  that  th 
volume  and  bulk  modulus  of  an  alloy  are  independent  of  the  chemn  . 
order.  This  assumption  allows  a  simple  calculation  of  the  elastic  part 
of  the  energy.  G,  dependent  solely  upon  the  concentration  and  the 
volume-independent  chemical  part,  e,  dependent  upon  the 
configuration.  Sluiter  et  al  [13]  proposed  a  slightly  more  realistic 
model  that  is  rigorously  correct  if  the  atoms  in  the  alloy  were  jellium 
soft  spheres.  The  model  uses  a  relaxation  energy  that  is  quadratic  in 
the  volume.  This  energy  is  subtracted  from  the  formation  energy 
leaving  the  unrelaxed  energy.  This  unrelaxed  energy  is  used  in 
calculation  of  the  configurational  free  energy  using  the  Cluster- 
Variation  Method  (CVM),  In  contrast  to  the  two  preceding 
techniques,  which  account  solely  for  volumetric  relaxation,  Zunger  et 
al  [24]  included  some  structural-relaxation  modes  in  the  total-energy 
calculations,  thereby  including  more  realistic  relaxation  energies  in 
the  configurational  free  energy. 

Except  in  the  thermodynamic  limit  of  the  expansion,  the 
interaction  energies  obtained  from  the  inversion  method  are  not 
necessarily  invariant  irrespective  of  the  choice  of  th°  ordered 
structures  used  in  the  expansion.  Nor  is  there  an  a  priori  knowledge 
of  the  range  (both  length  and  complexity)  beyond  which  interactions 
become  insignificant.  These  are  the  two  principal  aspects  of 
convergence  of  the  cluster  expansion  and  as  such  are  examined 
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closely  in  this  (and  other)  first-principles  phase-stability  studies. 
However,  for  a  well  converged  cluster  expansion  the  interaction 
energies  are  sufficiently  precise  to  reproduce  formation  energies  within 
5%  of  the  total-energy  calculations. 

At  finite  temperatures,  disorder  and,  commensurately, 
configurational  entropy  become  significant.  The  Ising  model  lends 
itself  to  two  solution  methods:  Monte  Carlo  techniques  [25-28]  and 
the  CVM.  Several  studies  have  gauged  the  accuracy  of  the  CVM  by 
comparison  with  Monte  Carlo  calculations  (see.  for  example,  [28]). 

The  CVM  was  proposed  by  Kikuchi  [29],  simplified  by  Barker  [30],  and 
reformulated  by  Morita  [31]  in  the  context  of  a  variational  statement 
of  statistical  mechanics.  Sanchez  and  de  Fontaine  [22,  32]  developed 
the  method  into  a  formulation  suitable  for  the  study  of  alloys.  The 
method,  which  reduces  to  the  well  known  Bragg-Williams  or  Bethe 
approximations  in  a  single-site  or  pair  approximation,  prescribes  a 
series  of  increasingly  higher  order  approximations  of  the  free-energy 
functional  to  an  arbitrary  level.  The  functional  may  then  be 
minimized  with  respect  to  the  configurational  variables,  i.e.,  the 
cluster  probabilities  or,  alternatively,  the  correlation  functions.  This 
minimization  of  the  configurational  free  energy  requires  calculation  of 
the  values  of  the  correlation  function  of  all  sub-clusters  of  the 
maximum  cluster.  A  small  increase  in  the  maximum  cluster  size 
from,  say,  the  joint  tetrahedron-octahedron  cluster  in  the  fee  lattice 
results  in  a  large  increase  in  the  number  of  sub-clusters,  thus 
increasing  significantly  the  number  of  operations.  As  with  the  cluster 
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expansion  of  the  formation  energy,  the  convergence  of  this  statistical 
approximation  must  be  assured.  In  order  to  include  longer-range 
interactions  in  the  Ising  Hamiltonian  without  increasing  the 
maximum  cluster  size,  many  studies  (for  example  (6,  10.  12,  16,  and 
33])  have  employed  the  so-called  renormalization  method  in  which 
long-range  interactions  are  included  using  a  Bragg-Williams  type 
approximation. 

Vibrational  free  energy  is  a  critical  driving  force  behind  some 
phase  transformations.  As  noted  in  Reference  15  the  leading  effect  of 
vibrational  entropy  is  responsible  for  the  curvature  of  the  fcc-bcc 
phase  boundary  in  the  Cu-Zn  phase  diagram.  However,  ab  initio 
treatments  of  vibrational  effects  are  rarely  included  in  first-principles 
phase-stability  studies.  Instead,  either  alloy  systems  are  chosen  so 
that  the  vibrational  free  energy  is  expected  to  be  insignificant  (e.g., 
all  structures  having  same  parent  Bravais  lattice),  or  atomic 
vibrations  are  included  empirically.  Since  ab  initio  calculations  of 
phonon  spectra  of  compounds  are  common  ]e.g.  34],  a  state-of-the-art 
first-principles  study  should  incorporate  vibrational  effects  without 
resort  to  measurements.  Moruzzi  et  al  [35]  calculated  Debye 
temperatures  of  several  cubic  transition  metals  using  a  very  simple 
model  with  good  results.  Their  model,  based  on  an  effective  speed  of 
sound  derived  from  the  bulk  modulus,  incorporates  quasi-harmonic 
dependence  of  the  vibrational  free  energy  through  the  Debye- 
Gruneisen  approximation.  This  technique  is  well  suited  to  the  CVM 
and  is  employed  in  this  study. 
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The  first  chapter  of  this  thesis  summarizes  the  reported 
thermodynamic  measurements  of  the  three  binary  systems  (Zr-Nb,  Nb- 
Ru.  and  Zr-Ru),  and  the  results  of  the  total-energy  calculations  using 
the  Linear  Muffin-Tin-Orbital  Method  in  the  Atomic-Sphere 
Approximation  (LMTO-ASA)  are  presented.  Reported  experimental 
data  for  the  systems  range  in  abundance  for  the  three  systems.  For 
Zr-Nb,  phase  diagram  features  and  calorimetry  data  have  been  well 
documented.  The  phase  diagram  for  Zr-Ru  has  been  established,  but 
few  measurements  have  been  reported  for  the  equilibrium  phases 
other  than  solubility  limits.  The  phase  diagram  for  Nb-Ru  is  by  no 
means  determinate.  Studies  of  the  crystallography  have  pointed  to 
the  phase  diagram  compiled  by  Massalsld  [36]  (see  Figure  1.5),  but 
results  from  the  present  study  and  those  of  Raub  and  Fritzsche  [37] 
and  Das  et  al  [38]  show  that  diagram  to  be  incomplete. 

Formation  energies  are  extracted  from  the  total-energy 
calculations  and  are  approximated  using  a  volume-dependent  Morse 
potential.  Several  ground-state  properties  such  as  bulk  moduli  and 
lattice  constants  are  then  available  from  these  fit  potentials  which 
are  compared  with  reported  measurements.  The,  as  yet,  undetermined 
crystal  structure  of  the  NbRu3  phase  is  posited  to  be  hexagonal 
DO i9.,  and  a  "newly  discovered"  stable  ground-state  phase,  NbsRu 
with  bcc  DO3  structure  is  predicted. 

Chapter  2  contains  a  concise  description  of  the  cluster 
expansion  method  and  an  in-depth  analysis  of  its  implementation  in 
this  study.  Several  optimization  schemes  are  examined,  with  no  one 
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method  emerging  as  definitive  --  though  they  prove  to  be  invaluable. 
The  convergence  of  the  cluster  expansion  is  gauged  according  to  how 
precisely  the  formation  energy  of  an  ordered  configuration  can  be 
estimated  using  the  derived  cluster  interactions.  Through  this  energy 
expansion  the  ground-state  crystal  structures  of  the  phases  NbaRu 
and  NbRu3  are  determined  to  be  the  most  stable. 

Finite-temperature  free-energy  calculations  are  presented  in 
Chapter  3.  The  CVM  expression  for  the  configurational  free  energy  is 
set  forth,  and  the  solution  method  developed  by  Sanchez  and  de 
Fontaine  [32J  is  outlined.  Included  is  the  treatment  of  vibrational 
free  energy  after  Moruzzi  [351.  This  allows  calculation  of  thermal 
expansion,  and  the  coefficients  of  thermal  expansion  of  the  pure 
elements  as  a  function  of  temperature  are  calculated  and  compared 
with  measurements.  The  calculated  miscibility  gap  in  Zr-Nb  is  in 
excellent  agreement  with  measurements. 


CHAPTER  1.  REPORTED  MEASUREMENTS  AND  TOTAL  ENERGY 
CALCULATIONS 


In  this  chapter  a  review  of  reported  measurements  of  thermochemical 
properties  of  the  three  binary  systems  Zirconium-Niobium,  Niobium- 
Ruthenium,  and  Zirconium-Ruthenium  is  provided.  In  Section  1.5 
the  results  of  the  total-energy  calculations  are  presented,  and  some  of 
the  calculated  properties  are  compared  with  measurements. 


l.l  Elemental  properties 

Table  1.1  includes  several  basic  properties  of  the  constituent 
elements  of  the  system.  Notably,  the  bulk  modulus  and  density 
increase  with  atomic  number,  but  the  melting  temperature  does  not. 
The  hep  -»  bcc  transition  of  zirconium  at  866  °C  is  driven  primarily 
by  vibration  modes  (39). 


1.2  Zirconium-Niobium 

Zirconium-Niobium  is  a  simple  system  with  small,  positive 
formation  enthalpies.  Figure  1.1  shows  the  phase  diagram  as 
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Table  1.1.  Measured  properties  of  Zr,  Nb,  and  Ru. 


Temperature 


compiled  by  Gullermet  {411.  The  hlgh-temperature  miscibility  gap  in 
the  bcc  phase  field  has  a  critical  temperature  reportedly  at  988  °C  at 
an  approximate  concentration  of  60.6  %  Nb  [421.  The  monotectoid 
equilibrium  temperature  Is  reported  as  610-620  °C  with  (Zr) 
concentration  between  0.6  and  1  weight  percent  Zr,  Zr-rich  bcc 
concentration  between  15  and  20.3  weight  percent  Zr,  and  Nb-rich 
bcc  concentration  of  85  -  93  weight  percent  Zr  [42J. 


1.3  Zirconium-Ruthenium 

Zirconium-Ruthenium  Is  also  a  relatively  simple  system.  Figure  1.2 
shows  the  assessed  phase  diagram  of  Raub  and  Rdschel  [43J.  The 
salient  feature  of  the  system  is  the  compound  ZrRu  with  B2  (CsCl) 
structure.  The  ordered  phase  is  stable  until  melting  at  approximately 
2100  °C.  A  hlgh-temperature  (above  1300  °C)  Laves  phase  (C14)  at 
ZrRu2  Is  destabilized  at  temperatures  below  1285  °C. 


Temperature 


u 

a 
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Weight  Percent  Zirconium 


Figure  1.2.  Zr-Ru  phase  diagram  from  Reference  43. 
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1.4  Niobium-Ruthenium 

The  niobium-ruthenium  phase  diagram  is,  as  yet, 

Incomplete.  In  fact,  some  controversy  exists  among  authors 
concerning  the  pure  ruthenium  phase.  Jaeger  and  Rosenboom  [44] 
reported  studies  suggesting  five  allotropic  transformations  in 
ruthenium.  However,  Hall  and  Crangle  [45]  performed  high- 
temperature  x-ray  measurements  showing  a  continuous  phase  with 
the  Mg  structure.  Their  conclusion  is  that  the  anamolous  calorimetry 
data  are  due  to  non-structural  second-order  transitions  of  unspecified 
nature.  Nevertheless.  Panteleemenov  et  al  [46]  proposed  a  phase 
diagram  for  Nb-Ru  assuming  the  existence  of  the  high-temperature 
polymorph  8. 

Greenfield  and  Beck  [47J  made  diffraction  measurements  of 
several  alloys  in  the  system.  They  postulated  a  B2-type  ordering  of 
the  near  equi-atomic  alloys  based  on  the  report  of  Dwight  [48],  but  the 
scattering  power  of  Nb  and  Ru  are  nearly  indistinguishable  for 
Chromium-K  X-rays,  and  no  superlattice  peaks  could  be 
distinguished.  However,  vanadium  and  tantalum  which  are  in  the 
same  group  as  zirconium,  form  B2  phases  with  ruthenium. 

Raub  and  Fritzsche  [37]  and  Hurley  and  Brophy  [49] 
reported  a  tetragonal  distortion  at  a  concentration  of  approximately 
40%  Ru.  Raub  and  Fritzsche  postulated  narrow  two-phase  regions  of 
bcc  +  bet  near  40  atomic-percent  Ru  and  bet  +  rhombohedral  near  46 
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atomic-percent  Ru.  Hurley  and  Brophy  mistakenly  concluded  that  a 
single  (Nb)  phase  field  spans  concentrations  up  to  55  atomic -percent 
Ru.  and  that  the  equi-atomic  alloy  Is  a  "distorted  bcc  phase."  Das  et 
al  [38]  also  reported  the  tetragonal  phase  at  concentrations  up  to  46 
%  Ru  where  an  orthorhombic  distortion  occurs.  Figure  1.3  shows  the 
lattice-parameter  measurements  of  Das  et  al  and  Raub  and  Fritzsche 
of  Nb-Ru  alloys  as  a  function  of  Ru  concentration.  It  is  clear  that  both 
studies  are  in  close  agreement  concerning  the  tetragonal  distortion, 
but  they  differ  as  to  the  lattice  type  of  the  equi-atomic  alloy.  Das  et  al 
determined  the  transition  temperature  of  the  tetragonal  phase  to  the 
parent  cubic  phase  at  three  concentrations  (45.8%.  51.1%.  and 
55.8%  Ru)  using  resistance  measurements.  Based  on  their 
measurements  and  the  reports  of  Raub  and  Fritzsche  and  of 
Greenfield  and  Beck  they  posited  a  low-temperature  ordered  face- 
centered  orthorhombic  phase,  a  higher-temperature  face-centered 
tetragonal  phase  with  Llo  ordering,  and  the  high-temperature  parent 
cubic  phase  with  B2  type  ordering.  However,  no  direct  evidence  of 
ordering  is  reported  due  to  the  difficulty  of  detecting  superlattice  lines 
for  the  system.  Tsukamoto  et  al  [50]  also  measured  the  temperature 
of  the  tetragonal  to  cubic  transition.  Figure  1.4  shows  the  phase 
boundary  between  the  cubic  and  tetragonal  phases  for  Nbso+xRuso-x 
proposed  by  Tsukamoto  et  al 

Hurley  and  Brophy  reported  a  hexagonal  compound  of 
unknown  structure  at  NbRu3.  Popova  [51]  synthesized  a  high- 
pressure,  high-temperature  compound  with  the  LI  2  (CU3AU)  structure 


Figure  1.3.  Measured  lattice  parameters  of  near  equi-atomic 
alloys  in  Nb-Ru.  From  Das  et  al  [38]  and  Raub  and  Fritzsche  [37]. 


Figure  1 .4.  Transition  temperature  of  L10  phase,  Nb^Ru^,. 
as  a  function  of  x.  (From  Tsukamoto  et  al  (50]). 
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at  NbRu3  (which,  incidentally,  is  superconducting  at  15  K);  however 
no  details  were  included  in  the  abstract  (51).  Raub  and  Fritzsche 
detected  no  distinct  phase  between  the  orthorhombic  phase  and  the 
terminal  Ru  phase. 

Figure  1.5  shows  the  proposed  phase  diagram  compiled  by 
Massalski  (361.  This  phase  diagram  is  essentially  identical  to  that  of 
Hurley  and  Brophy  except  for  the  tetragonal  (NbRu')  to  cubic  (NbRu) 
transition  boundary  reported  by  Tsukamoto  and  the  B2  +  bcc  two- 
phase  region  reported  by  Raub  and  Fritzsche.  Massalski  also  labels 
the  NbRu3  phase  as  LI  2  following  Popova.  No  reference  to  the 
orthorhombic  phase  reported  by  Raub  and  Fritzsche  and  by  Das  et  al 
is  included. 


Temperature 


Weight  Percent  Ruthenium 

10  20  30  40  50  60  70  80  90  100 


Figure  1.5.  Proposed  Nb-Ru  phase  diagram  from  Reference  36. 
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The  Hamiltonian  for  a  system  with  many  interacting  electrons  may  be 
written 

H  =  T  +  U  +  V,  (1.1] 

where  T  is  the  kinetic  energy  operator.  U  is  the  Coulombic  electron- 
electron  potential,  and  V  is  the  external  potential  including  the 
electrostatic  interaction  with  fixed  nuclei.  The  Density  Functional 
Theory  (DFT)  [1]  states  that  V  is  a  unique  functional  of  the  electron 
density  and  that  the  minimum  energy  is  reached  for  the  correct 
electron  density.  The  many-electron  problem  is  intractable,  and  thus 
the  functional  may  not  be  ascertained.  Kohn  and  Sham  [2]  applied 
the  DFT  formalism  to  the  one-electron  problem  by  including  a  suitable 
functional  for  the  external  potential.  This  functional  requires  the 
calculation  of  the  one-electron  kinetic  energy  and  an  exchange- 
correlation  energy.  The  Local-Density  Approximation  (LDA)  is  made 
by  adopting  an  analytical  form  of  the  exchange-correlation  functional 
which  depends  on  the  electron  density  and  is  exact  for  a 
homogeneous  electron  gas.  The  calculation  of  the  one-electron  band 
structure  is  then  required  to  obtain  the  kinetic  energy.  For 
computational  expediency  linear  band-structure  methods  are 
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required.  Such  linear  techniques  as  FLAPW  (full-potential  augmented 
plane-wave  [52]).  ASW  (augmented  spherical  wave  [531),  and  LMTO 
(linearized  muffin-tin  orbital  (54,  20))  are  well-developed  and  in 
general  use.  In  the  LMTO  method,  energy-independent  basis 
functions  are  constructed  as  the  solutions  of  the  one-electron 
SchrOdinger  equation  for  the  atomic-like  muffin-tin  potential.  The 
solution  method  has  leading  error  terms  of  second  order  in  energy.  In 
this  study,  the  Atomic-Sphere  Approximation  (ASA),  which  restricts 
the  muffin-tin  potentials  to  spherical  symmetry,  was  included.  In  the 
ASA.  only  high-symmetiy  close-packed  structures  may  be  treated 
accurately. 

In  this  study,  spin-orbit  interactions  were  included  and  the 
exchange-correlation  potential  of  Hedin  and  Lunqvist  (55)  was  used. 
The  total  energies,  Etotal,  of  each  compound  were  calculated  at  12 
unit-cell  volumes  close  to  the  equilibrium  volume.  The  crystal 
structures  used  are  represented  in  Figures  1.6  a.  b.  and  c.  The 
figures  show  the  basal-plane  projections  of  one  or  more  unit  cells  of 
each  of  the  structures  in  the  bcc  and  hep  lattices.  Filled  and  empty 
circles  represent  atoms  of  type  a  or  b,  respectively.  The  smaller 
circles  denote  atoms  with  positions  displaced  one  half  a  unit  cell 
vector  normal  to  the  basal  plane.  Half-filled  circles  denote 
alternations  in  the  stacking  sequence  of  atoms  in  the  direction  normal 
to  the  basal  plane. 


T1 


T2 


Figure  1 .6  a.  Basal-plane  projections  of  bcc  based  ordered  structures  used 
in  this  study. 
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LiRh 


do19 


Figure  1.6  b.  Basal-plane  projections  of  hep  based  crystal 
structures  used  in  this  study.  The  smaller  circles  indicate 
lattice  positions  c/2  above  basal  plane. 
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The  cohesive  energy  per  atom,  Ecoh  .  of  each  compound 
AnBm  was  determined  at  each  unit  cell  volume  from: 


E^lAnBJ  = 


—total ,  .  _  ,  —atom  —atom 

E  (AnBm)  -  n-  Ea  -  m  EB 
m  +  n 


11.21 


where  Eatom  is  the  energy  of  the  free  atom  calculated  using  the 
Hartree  method.  The  formation  energy,  E,  is  then  calculated 
according  to  the  definition: 


coh 

E(AnBm)  =  E  (AnBm)  - 


_  ref  _  ref 

n-EA  +  m-Ee 
m  +  n 


11.3] 


where  Eref  is  the  reference  energy  of  each  pure  solid.  In  this  study 
the  cohesive  energy  of  each  element  in  its  stablest  configuration  Is 
used  as  the  reference  energy,  i.e.,  =  E^,  E^  =  E^,  and 

Ej^  =  E^ ,  unless  otherwise  noted.  Equation  1.3  may  be  generalized 
to  any  configurational  property  F(AnBm),  thus  defining  an  excess 
configurational  property: 


F^AnBm)  =  F(AnBm)  - 


—  ref  —  ref 

n-FA  +  m-Fe 
m  +  n 


11.4] 


A  Morse-type  potential  function  of  average  atomic  radius 
(Wigner-Seitz  radius),  r: 
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E(r)  =  A  -  2  D  exp(-\(r-ro))  +  D  exp(-2A.(  r-  ro))  ,  {1.51 

was  then  fit  to  the  formation  energies  for  each  compound.  Note  that 
the  minimum  energy  is  given  by, 

E  =  A  -  D  .  forr  =  ro.  [1.6] 

Also  note  that  if  the  cohesive  energy  of  a  solid  is  fit  with  a  Morse 
potential,  the  energy  approaches  that  of  the  free  atom(s)  in  the  large 
limit  of  r: 


E  — » A  ,  as  r  .  [1.7] 

Thus  the  parameter  A  should  be  a  linear  function  of  the  composition, 
and  there  are  only  three  degrees  of  freedom,  ro.  A.,  and  D,  in  the  fit. 
However,  due  to  (i)  the  artificial  shape  of  the  Morse  potential,  (ii)  the 
inaccuracy  inherent  in  the  solution  method  of  the  Hartree  method  for 
calculating  the  energy  of  the  free  atoms,  and  (ill)  the  fact  that  the 
Morse  function  is  fit  to  data  within  a  relatively  narrow  domain  of  r,  the 
fit  values  of  A  do  not  agree  precisely  with  the  calculated  free-atom 
energies  (errors  are  typically  10-20  mRy/atom).  This,  however,  is  not 
a  significant  problem,  as  the  fit  of  the  potential  within  the  domain  of 
interest  --  atomic  radii  between  those  of  the  constituent  atoms  --  has 
discrepancies  on  the  order  of  tenths  of  milli-Rydbergs.  Figure  1.7 
shows  the  Morse  fit  to  the  calculated  formation  energies  of  NbRu3 
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with  a  DO  19  structure.  The  Inset  figure  shows  the  overall  shape  of  the 
Morse  potential  and  indicates  the  scale  of  the  formation  energy  with 
respect  to  the  cohesive  energy. 

The  calculated  Morse  parameters  and  minimum  energies  for 

each  compound  studied  are  shown  in  Table  1.2  a-c.  The  average 
atomic  volume  is  Do  =  fn  ro  .  The  experimentally  determined 

o 

equilibrium  atomic  volumes  of  the  pure  elements  compare  with  those 
calculated  as  follows: 

hcp-Zr  experimental  *  157.00  a.u.3  calculated  -  166.82  a.u.3 

hcp-Ru  experimental  -  91.85  a.u.3  calculated  -  99.78  a.u.3 

bcc-Nb  experimental  -  133.18  a.u.3  calculated  -  121.26  a.u.3 

Typically.  LDA  calculations  predict  lower  equilibrium  volumes  than 
measurements,  and  this  effect  is  seen  for  niobium.  However  the 
calculations  for  hexagonal  structures  were  confined  to  lattices  with 
ideal  close-packing,  i.e.,  c/a  =  (8/3)1/2  »  1.63,  whereas  the  measured 
c/a  ratios  for  Zr  and  Ru  from  Table  1.1  are  1.59  and  1.58  respectively. 
This  artificially  high  value  for  c/a  accounts  partially  for  the 
overestimation  of  the  atomic  volume  in  the  hexagonal  structures.  Due 
to  the  constraints  of  the  ASA  the  total-energy  calculations  were 
generally  confined  to  close-packed  structures  in  this  study  —  though 
some  results  for  distorted  lattices  were  included  and  are  discussed 
later  in  this  section.  Lu  et  al  [56]  used  the  (FLAPW)  method  to 


Table  1.2  a.  The  Morse  parameters  fit  to  excess  structural  energy  for  each 
compound  In  Nb-Ru  studied. 
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|  Compound/ 
Structure 

ro 

(a.u.) 

\ 

(a.u.-1! 

D 

[Ry.l 

A 

(Ry.l 

Emin 

(Ry.l 

Nb/hcp 

3.1878 

1.0423 

0.71738 

0.74348 

0.026101 

NbsRu/DOig 

3.1046 

1.109 

0.7056 

0.71511 

0.009519 

NbRu/Pmma 

3.0267 

1.1887 

0.69147 

0.68674 

-0.00473 

NbRus/DOjg 

2.9494 

1.2793 

0.67003 

0.65837 

-0.01166 

Ru/hcp 

2.8773 

1.3938 

0.63 

0.63 

0 

NbRu/LiRh 

3.0228 

1.178 

0.68868 

0.68674 

-0.00193 

NbRu/Pmmn 

3.0273 

1.1825 

0.68685 

0.68674 

-0.0001 

NbRu/AB5 

3.026 

1.186 

0.68989 

0.68674 

-0.00314 

Nb/fcc 

3.1912 

1.0451 

0.71835 

0.74348 

0.025128 

Nb3Ru/Ll2 

3.1071 

1.1105 

0.70641 

0.71511 

0.008705 

NbRu/Llo 

3.0265 

1.1905 

0.69351 

0.68674 

-0.00676 

NbRu3/Ll2 

2.9492 

1.2808 

0.66935 

0.65837 

-0.01098 

Ru/fcc 

2.8813 

1.3936 

0.62399 

0.63 

0.006015 

Nb3Ru/D022 

3.1066 

1.1106 

0.70559 

0.71511 

0.009521 

NbR.U3/DC>22 

2.9519 

1.2824 

0.6664 

0.65837 

-0.00802 

NbRu/NbP 

3.0263 

1.1938 

0.69207 

0.68674 

-0.00532 

NbRu/Tl 

3.0238 

1.1966 

0.69455 

0.68674 

-0.00781 

Nb/bcc 

3.1678 

1.0438 

0.74348 

0.74348 

0 

Nb3Ru/D03 

3.0887 

1.1153 

0.73233 

0.71511 

-0.01722 

NbRu/B2 

3.0313 

1.1966 

0.69953 

0.68674 

-0.01279 

NbRu/B32 

3.0261 

1.1845 

0.69221 

0.68674 

-0.00546 

NbRu3/D03 

2.9698 

1.2769 

0.64682 

0.65837 

0.011555 

Ru/bcc 

2.9095 

1.3885 

0.58893 

0.63 

0.041073 

NbRu/Tl 

3.025 

1.187 

0.69669 

0.68674 

-0.00995 

NbRu/T2 

3.0329 

1.1987 

0.69237 

0.68674 

-0.00563 

NbRu/T3 

3.0329 

1.1987 

0.69235 

0.68674 

-0.00561 
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Table  1.2  a  continued 


r£) 

D 

A 

Emin 

NbRu/T4 

3.0334 

1.1856 

0.69088 

0.68674 

-0.00414 

Nb3Ru/T5 

3.0968 

1.1148 

0.73046 

0.71511 

-0.01535 

NbRu3/T5 

2.9728 

1.2805 

0.64568 

0.65837 

0.01269 

Nb/bcc  * 

3.1725 

1.0388 

0.73884 

0.74348 

0.004647 

Nb3Ru/D03* 

3.094 

1.1209 

0.72558 

0.71511 

-0.01047 

NbRu/B2  * 

3.0278 

1.2062 

0.70209 

0.68674 

-0.01534 

NbRu/B32  • 

3.026 

1.1934 

0.69026 

0.68674 

-0.00351 

NbRu3/D03* 

2.9649 

1.2825 

0.65243 

0.65837 

0.005945 

Ru/bcc  • 

2.9005 

1.3878 

0.60048 

0.63 

0.029522 

NbRu/Tl  * 

3.0255 

1.1909 

0.69489 

0.68674 

-0.00814 

NbRu/T2  • 

3.0324 

1.1939 

0.6927 

0.68674 

-0.00596 

NbRu/T3  • 

3.0333 

1.1898 

0.68803 

0.68674 

-0.00128 

NbRu/T4  * 

3.0368 

1.1898 

0.68718 

0.68674 

-0.00043 

Nb3Ru/T5  * 

3.1005 

1.1197 

0.72544 

0.71511 

-0.01032 

NbRu3/T5  • 

2.9701 

1.2814 

0.64756 

0.65837 

0.010814 

•  denotes  distorted  bcc  lattice  with  c/a=  1.122. 


Table  1.2  b.  The  Morse  parameters  fit  to  excess  structural  energy  for  each 
compound  in  Zr-Ru  studied. 
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Compound/ 

Structure 

r0 

la.u.l 

X 

(a.u.-1! 

D 

[Ry.] 

A 

[Ry.) 

Emin 

[Ry.l 

Zr/bcc 

3.3994 

0.84306 

0.68675 

0.68923 

0.002483 

Zt3Ru/DC>3 

3.2437 

0.963 

0.69483 

0.67442 

-0.02041 

ZrRu/B2 

3.1221 

1.0892 

0.71137 

0.65962 

-0.05175 

ZrRu/B32 

3.1117 

1.0636 

0.6806 

0.65962 

-0.02099 

ZrRu3/D03 

3.014 

1.1992 

0.64553 

0.64481 

-0.00072 

Ru/bcc 

2.9095 

1.3885 

0.58893 

0.63 

0.041073 

ZrRu/Tl 

3.1125 

1.0593 

0.67913 

0.65962 

-0.01951 

ZrRu3/T2 

3.1206 

1.0617 

0.67744 

0.65962 

-0.01/82 

ZrRu/T3 

3.1191 

1.0577 

0.67142 

0.65962 

-0.0118 

ZrRu/T4 

3.1203 

1.0458 

0.66821 

0.65962 

-0.0086 

Zr3Ru/T5 

3.2536 

0.94904 

0.68999 

0.67442 

-0.01557 

ZrRu3/T5 

3.0155 

1.2101 

0.64394 

0.64481 

0.000865 

Zr/hcp 

3.4154 

0.85552 

0.68923 

0.68923 

0 

Zr3Ru/DOi9 

3.2654 

0.94832 

0.67978 

0.67442 

-0.00535 

ZrRu/Pmma 

3.127 

1.0589 

0.67989 

0.65962 

-0.02027 

&RU3/DO19 

3.004 

1.1976 

0.65615 

0.64481 

-0.01134 

Ru/hcp 

2.8773 

1.3938 

0.63 

0.63 

0 

ZrRu/LiRh 

3.1248 

1.0405 

0.66751 

0.65962 

-0.00789 

ZrRu/Pmmn 

3.1264 

1.0475 

0.66679 

0.65962 

-0.00718 

ZrRu/Hl 

3.1252 

1.0645 

0.67235 

0.65962 

-0.01273 

Table  1.2  c.  The  Morse  parameters  fit  to  excess  structural  energy  for  each 
compound  In  Zr-Nb  studied. 
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Compound/ 

Structure 

ro 

(a.u.J 

31 

lau.'1] 

D 

IRy.] 

A 

IRy.) 

Emin 

[Ry] 

Zr/bcc 

3.399 

.84306 

0.68675 

0.68923 

0.00252 

Zr3Nb/D03 

3.337 

.8886 

0.69913 

0.70279 

0.00366 

ZrNb/B2 

3.281 

.9366 

0.70905 

0.71635 

0.00730 

ZrNb/B32 

3.278 

.9397 

0.71346 

0.71635 

0.00289 

ZrNb3/D03 

3.222 

.9903 

0.72665 

0.72992 

0.00326 

Nb/bcc 

3.168 

1.0438 

0.74348 

0.74348 

0.00000 

ZrNb/Tl 

3.277 

.9386 

0.71182 

0.71635 

0.00453 

ZrNb/T2 

3.278 

.9366 

0.71047 

0.71635 

0.00588 

ZrNb/T3 

3.277 

.9367 

0.71104 

0.71635 

0.00532 

ZrNb/T4 

- 

- 

0.71013 

0.71635 

0.00622 

Zr/hcp 

3.415 

.8554 

0.68923 

0.68923 

0.0000 

Zr3Nb/DOi9 

3.360 

.8843 

0.69120 

0.70279 

0.0116 

ZrNb/Pmma 

3.302 

.9358 

0.69678 

0.71635 

0.0196 

ZrNb/Pmmn 

3.302 

.9418 

0.69705 

0.71635 

0.0193 

ZrNb3/DOi9 

3.246 

.9870 

0.70509 

0.72992 

0.0248 

Nb/hcp 

3.188 

1.0423 

0.71738 

0.74348 

0.0261 

ZrNb/Hl 

3.304 

.9381 

0.69965 

0.71635 

0.0167 

ZrNb/LIRh 

3.303 

.9438 

0.69718 

0.71635 

0.0192 

calculate  total  energies  for  hcp-Zr.  They  predicted  an  equilibrium 
volume  of  158.92  a.u.3  and  a  c/a  ratio  of  1.593. 


The  pressure  P  is  simply: 


and  the  bulk  modulus  B  is: 

B=-vi=v0-  1191 

or,  at  the  equilibrium  volume 

2 

B(r0)  =  r  11.10] 

6rcr0  6.7977  E-9  (Ry/a.u.3) 


The  calculated  bulk  moduli  of  each  compound  are  shown  in  Tables 
1.3  a-c.  The  only  measured  values  of  bulk  moduli  reported  for  these 
systems  are  for  the  pure  elements  (see  Table  1.1).  The  calculated  bulk 
modulus  of  Niobium  is  1.996  Mbar,  the  reported  value  1.7  Mbar.  For 
Ruthenium  the  calculated  bulk  modulus  is  3.319  Mbar,  the  measured 
value,  3.21  Mbar.  For  Zirconium  the  calculated  value  is  1.153  Mbar, 
the  reported  value  is  0.83  Mbar. 
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Table  1.3  a.  Calculated  bulk  moduli  for  each  compound  in  Nb-Ru  studied. 


Compound/ 

Structure 

bulk 

modulus 

JMbarl 

Nb/hcp 

1.908 

Nb3Ru/DOi9 

2.181 

NbRu/Pmma 

2.519 

NbRu3/DOi9 

2.896 

Ru/hcp 

3.319 

NbRu/LiRh 

2.467 

NbRu/Pmmn 

2.475 

NbRu/Hl 

2.503 

Nb/fcc 

1.919 

Nb3Ru/Ll2 

2.188 

NbRu/Llo 

2.535 

NbRu3/Ll2 

2.906 

Ru/fcc 

3.283 

Nb3Ru/D022 

2.186 

NbRu3/D022 

2.898 

NbRu/NbP 

2.543 

Nb/bcc 

1.996 

Nb3Ru/D03 

2.301 

NbRu/B2 

2.579 

NbRu/B32 

2.505 

Compound/ 

Structure 

bulk 

modulus 

[MbarJ 

NbRu3/DC>3 

2.771 

Ru/bcc 

3.046 

NbRu/Tl 

2.532 

NbRu/T2 

2.560 

NbRu/T3 

2.560 

NbRu/T4 

2.499 

Nb3Ru/T5 

2.288 

NbRu3/T5 

2.780 

Nb/bcc  • 

1.916 

Nb3Ru/DC>3* 

2.300 

NbRu/B2  * 

2.633 

NbRu/B32  * 

2.536 

NbRu3/D03* 

2.825 

Ru/bcc  * 

3.112 

NbRu/Tl  • 

2.542 

NbRu/T2  * 

2.541 

NbRu/T3  * 

2.506 

NbRu/T4  * 

2.500 

Nb3Ru/T5  * 

2.289 

NbRu3/T5  * 

2.794 

*  denotes  distorted  bcc  lattice  with  c/a=1.122. 
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Table  1.3  b.  Calculated  bulk  moduli  for  each  compound  in  Zr-Nb  studied. 


Compound/ 

bulk 

Structure 

modulus 

[Mbarl 

Zr/hcp 

1.153 

&3Nb/DO  19 

1.282 

ZrNb/Pmma 

1.454 

ZrNb3/DOi9 

1.660 

Nb/hcp 

1.908 

ZrNb/LiRh 

1.463 

ZrNb/Pmmn 

1.460 

ZrNb/Hl 

1.454 

Zr/bcc 

1.120 

Zr3Nb/D(>3 

1.291 

ZrNb/B2 

1.479 

ZrNb/B32 

1.500 

ZrNb3/D03 

1.726 

Nb/bcc 

1.996 

ZrNb/Tl 

1.493 

ZrNb/T2 

1.484 

ZrNb/T3 

1.486 

35 


Table  1.3  c.  Calculated  bulk  moduli  for  each  compound  In  Zr-Ru  studied. 


Compound/ 

Structure 

bulk  modulus 

[Mbar) 

Zr/hcp 

1.153 

ZraRu/DOig 

1.461 

ZrRu/Pmma 

1.902 

ZrRu3/DOi9 

2.445 

Ru/hcp 

3.319 

ZrRu/LiRh 

1.875 

ZrRu/Pmmn 

1.826 

ZrRu/Hl 

1.854 

Zr/bcc 

1.120 

Zr3Ru/D03 

1.552 

ZrRu/B2 

2.102 

ZrRu/B32 

1.933 

ZrRu3/D03 

2.406 

Ru/bcc 

3.046 

ZrRu/Tl 

1.914 

ZrRu/T2 

1.912 

ZrRu/T3 

1.882 

ZrRu/T4 

1.830 

Zr3Ru/T5 

1.492 

ZrRu3/T5 

2.435 
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Zr-Ru 


The  formation  energies  of  the  compounds  in  Zr-Ru  are 
displayed  graphically  in  Figure  1.8.  The  only  stable  intermediate 
compound  ZrRu  has  a  B2  (CsCl)  structure.  The  calculated  formation 
enthalpy  is  51.75  milli-Rydbergs  per  atom.  The  measured  value  is 
52.29  ±  2.59  mRy/atom  [57].  The  ground-state  line  extends  from  the 
pure  element  points  --  which  by  definition  are  zero  —  through  the 
compounds  with  the  lowest  energies  at  each  concentration.  In  Figure 
1.8  no  compound  has  an  energy  which  is  lower  than  the  linear 
combination  of  the  B2  compound  with  one  of  the  pure  elements.  Of 
course,  only  a  small  number  (11)  of  total  energy  calculations  have 
been  performed  for  this  system,  so  the  question  as  to  whether  an 
unexamined  compound  might  be  stable  looms  large.  This  question  is 
addressed  in  Chapter  2. 

Pearson  [58]  reported  the  average  atomic  volume  for  ZrRu 
(B2)  as  116.15  a.u.3,  the  computed  equilibrium  volume  is  127.48 
a.u.3.  Mehl  et  al  [59]  calculated  the  total  energy  versus  volume  for 
ZrRu  using  the  FLAFW  method  and  predicted  an  equilibrium  volume 
of  112.65  a.u.3. 

Figures  1.9  a  and  b  show  the  calculated  atomic  volumes  and 
bulk  moduli  of  each  compound  in  Zr-Ru.  The  excess  atomic  volume 
calculated  according  to  Equation  1.4  is  also  shown  in  Figure  1.9  a, 
and  the  dashed  lines  indicate  a  linear  rule  of  mixtures  for  each  lattice 
type.  The  average  volume  obeys  the  linear  rule  with  small  negative 


Ru 


Concentration  of  Zr 


Figure  1.8.  Calculated  formation  energies  for  selected  compounds 
in  Zr-Ru.  The  reported  value  for  the  B2  compound,  52.29  mRy/atom 
is  from  [57]. 


Excess  atomic  volume 
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deviations  (approximately  6%).  The  bulk  moduli  of  the  intermediate 
compounds  deviate  more  drastically  (18%),  and,  as  seen  from  Table 
1.3  c,  the  bulk  modulus  exhibits  significant  dependence  upon  the  type 
of  ordering.  The  calculated  bulk  modulus  of  ZrRu  (B2)  is  2. 102  Mbar, 
and  Mehl  et  al  (59]  predicted  2.25  Mbar. 


Nb-Ru 


The  ground-state  diagram  for  Nb-Ru  is  shown  in  Figure  1.10. 
The  ground-state  line  passes  through  Nb3Ru  with  the  DO3  structure. 
No  such  compound  has  been  reported.  Possibly  due  to  the  small 
scattering  difference  between  Nb  and  Ru.  it  is  difficult  to  discern  the 
ordered  bcc  structure  from  diffraction  measurements. 

A  stable  Llo  phase  with  c/a  =  1.122  is  seen  in  the  ground- 
state  diagram,  and  the  average  of  the  measured  c/a  and  c/a  ratios  of 
the  orthorhombic  phase  at  that  composition  is  1.1 165  110].  This  is 
remarkable  since  the  LMTO-ASA  typically  loses  accuracy  with 
deviations  from  close-packing.  Part  of  the  explanation  lies  in  the  fact 
that  at  c/a  =  1  the  lattice  is  bcc  and  at  c/a  =  V2=  1.414  the  lattice  is  fee 
(see  Figure  1.11),  thus  bounding  the  aspect  ratio  of  the  Llo  structure 
with  close-packed  lattices.  Further  remarks  on  the  stabilization  of  the 
Llo  phase  are  found  in  Section  2.6. 


Formation  energy  [mRy/atom] 


Figure  1.11.  Change  of  crystal  symmetry  due  to  tetragonal 
distortions  of  B2  structure. 
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The  stable  structure  of  the  NbRu3  compound  Is  not  obvious 
from  the  ground-state  diagram.  Table  1.2  a  shows  the  hexagonal 
DO  19  structure  to  have  a  slightly  lower  formation  energy  than  the  LI2 
structure.  Notably,  the  constraint  of  close-packing  in  the  total-energy 
calculations  yields  a  less  negative  energy  than  would  be  found  upon 
relaxation  of  the  c/a  ratio,  thus  the  DO  19  structure  is  certainly  the 
more  stable. 

The  atomic  volumes  and  bulk  moduli  of  the  fee-,  hep-,  and 
bcc-lattice  structures  in  Nb-Ru  are  shown  in  Figure  1.12.  As  with  Zr- 
Ru,  the  deviations  from  an  ideal  rule  of  mixtures  of  the  atomic 
volumes  is  negative  for  each  lattice  type.  That  is.  a  straight  line 
between  the  volumes  of  the  pure  elements  (indicated  by  the  dashed 
lines  in  Figure  1.12)  of  a  single  lattice  type  lies  above  the  volumes  of 
all  the  intermediate  compounds  of  the  same  lattice  type.  The 
deviations  of  the  volumes  and  of  the  bulk  moduli  from  the  straight 
line  mixture  rule  are  small  (less  than  3%). 

Zr-Nb 


The  ground-state  energies  of  the  Zr-Nb  system  are  shown  in 
Figure  1.13.  As  expected  from  the  well-documented  phase  diagram 
there  are  no  stable  intermediate  compounds.  Figure  1.14  shows  the 
calculated  atomic  volumes  and  bulk  moduli  of  each  compound.  As  in 
the  other  two  systems,  the  volumes  and  bulk  moduli  of  the 


46 


intermediate  compounds  have  small  negative  deviations  from  linear 
mixture  theory. 

The  total-energy  calculations  predict  ground-state  structures 
matching  those  reported  by  experimentalists  with  the  exception  of  the 
orthorhombic  phase  of  near  equi-atomic  Nb-Ru  alloys  which  is  not 
within  the  realm  of  the  atomic  sphere  approximation.  Additionally, 
the  ordering  of  the  compound  NbRu3  is  posited  to  be  DO  19,  and  an 
undetected  ground-state  phase  NbsRu  with  DO3  structure  is  also 
predicted.  The  calculated  atomic  volumes  of  the  elements  fall  within 
9%  of  the  measured  values.  The  calculated  atomic  volume  of  ZrRu  is 
within  10%  of  the  measured  value.  The  calculated  c/a  ratio  of  the 
NbRu  Llo  compound  is  within  0.5%  of  the  measured  value.  The 
calculated  bulk  moduli  of  the  Zr,  Nb,  and  Ru  are  within  18%  ••  and 
39%  respectively  of  the  reported  values,  and  the  predicted  valut  r 
ZrRu  is  within  7%  of  the  calculated  value  of  Mehl  et  al  [59]. 

The  total-energy  calculations  of  the  small  number  of 
structures  of  Table  1.1  are  not  sufficient  to  determine  the  complete 
ground-state  of  these  systems.  Due  to  the  fact  that  LDA  calculations 
for  structures  with  more  than  eight  atoms  per  primitive  unit  cell 
become  very  expensive,  only  a  small  select  set  of  structures  is 
practical.  The  cluster  expansion  described  in  Chapter  2  may  be  used 
to  estimate  the  formation  energy  of  any  ordered  structure  on  a  bcc. 
fee,  or  hep  lattices.  This  expansion  is  a  powerful  tool  which  lends 
confidence  to  the  predicted  ground-state  configurations  made  in  this 
chapter. 


CHAPTER  2.  CLUSTER  EXPANSION 


In  the  next  section,  highlights  of  the  general  formalism  of  the 
cluster  expansion  techniques  developed  by  Sanchez  et  al  [21,  22]  are 
reviewed.  Only  binary  systems  will  be  considered  here,  but  the  theory 
has  been  generalized  to  multi-component  systems  by  Sanchez  et  al 
[21].  The  main  result,  that  any  function  of  lattice  configuration 
(used  synonymously  in  this  chapter  with  configurational  property) 
can  be  described  by  a  cluster  expansion,  will  be  established,  as  will 
several  other  useful  notions. 

The  second  section  details  how  a  truncated  cluster 
expansion  can  be  inverted  to  derive  the  expansion  coefficients 
(effective  interactions)  from  a  set  of  known  configurational 
properties,  e.g.,  the  cohesive  energy.  Section  2.3  sets  forth  several 
schemes  for  including  structural  relaxation  in  the  cluster  expansion. 

A  truncated  cluster  expansion  is  not  unique,  and  the  final  three 
sections  of  this  chapter  describe  the  optimization  methods  used  to 
insure  adequate  convergence  of  the  expansion.  The  most  important 
method  for  testing  the  convergence  is  using  the  cluster  expansion  to 
predict  cohesive  energies  of  various  ordered  structures  for  which  the 
energies  have  been  previously  determined  using  the  local  density 
approximation  (LMTO)  calculations.  The  results  of  the  convergence 
tests  are  discussed  in  Section  2.6. 
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2J_^UgtgT.ggp{»iSl9n.&>rm.ftll8ni 

In  the  Islng  model  a  binary  crystal  with  N  lattice-sites  has  2N 
possible  configurations.  To  a  site  i  occupied  by  an  A-atom,  a  spin  si 
of  +1  is  assigned,  and  for  B-atoms  -1.  Each  of  the  2N  possible 
configurations  is  then  described  completely  by  the  N-vector  a  =  {  sj. 

S2.  S3 . sn)  .  A  set  of  orthonormal  characteristic  functions.  4>.  for 

N-dimensional  discrete  space  spanned  by  ct  is  given  by  (21J: 

=  n  s,  =  s„  s*  •  •  Sfc,  ,  [2.11 

lea 

where  a  =  { ii,  i2,  i3,  ...y  includes  every  cluster  of  n  sites  in  the 
crystal  lattice,  and  where  the  inner  product  is  defined  as: 

<Oa(a)-4>p(a )  >  =  <i>a(o)  Op(cy)  (2.2) 

2  ° 

Hence,  there  is  a  one-to-one  correspondence  between  the  set  of 
orthogonal  functions  O  and  the  set  of  all  clusters  a,  including  the 
empty  cluster  for  which.  C>o(o)  =  1. 

As  <I>  forms  a  complete  orthonormal  basis,  any  function  of 
configuration  (configurational  property).  F(c),  may  be  written  as  [21]: 
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F(o)=  IJ«<I>Jo)  .  [2.31 

a 

where  the  sum  extends  over  all  inequivalent  clusters  in  the  crystal, 
and  Ja  .  given  by: 

Ja=  <  F(a)-  <l>a(o)  > 

=  -VlF{o)  <Da(o)  [2.4J 

2  • 

are  the  projections  of  F(c)  on  the  orthogonal  cluster  basis,  <t>  —  often 
called  the  effective  interactions.  Equation  2.3  is  simplified  by  the 
requirement  that  each  Ja  be  the  same  for  all  clusters  that  are  related 
by  any  of  the  space-group  operations  of  the  crystal. 

The  correlation  functions,  2^,  are  defined  as  the  expectation 
values  of  the  characteristic  functions: 


%o.  =  <  Oja)  >  , 

=  <  SiS2--  -  Sna> 


[2.5] 


For  a  perfectly  ordered  crystal,  the  correlation  function  of  a  cluster 
type  is  an  average  of  the  spin  product  over  all  equivalent  clusters  in 
the  crystal: 
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4a  =  - ^ -  E  Spj  Spa  •  •  •  Sp„  (2.6] 

Za-N  W 

where  a  labels  each  type  of  Inequivalent  cluster  of  which  there  are 
z«-  N  in  the  crystal.  Combining  Equation  2.3  with  Equation  2.5 
leads  to  an  expression  for  the  expectation  value  of  any  configurational 
property  °F: 

«F=X<F(o)>=  i  Zji  Jn  4n  •  [2.7] 

N  n  -  0 

in  terms  of  the  configurational  variables.  4n.  and  the  interaction 
coefficients,  Jn.  The  technique  for  determing  the  interactions  is 
discussed  in  the  next  section. 


2.2  Obtaining  cluster  Interactions  from  configurational 
properties 

Since  the  sum  in  Equation  2.6  cannot  be  evaluated  exactly 
(there  are  on  the  order  of  1023  terms),  a  truncated  approximation 
must  be  used,  and  the  practical  application  of  the  representation  of 
the  expectation  value  of  a  configurational  function  depends  upon  the 
rapid  convergence  of  the  cluster  projections,  Jn.  The  key  step  is 
choosing,  a  priori,  a  trial  maximum  cluster  size  and  establishing  the 
convergence  of  the  truncated  series,  with  the  dual  task  of  determining 
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the  cluster  interactions.  The  former  is  discussed  in  section  2.6.  The 
latter  is  the  natural  first  step  and  is  accomplished  using  a  system  of 
equations  of  the  form  of  Equation  2.6  for  the  value  of  the  energy  of  an 
ordered  configuration  p: 

pE=Izkp$kJk.  [2.8] 

k*0 

where  k  denotes  the  cluster  type  and  is  the  correlation  function 
corresponding  to  cluster  k  in  the  ordered  structure  p  which  is 
evaluated  directly  from  the  geometry  of  the  configuration.  For 
convenience  Equation  2.8  is  rewritten: 

PE=  I  EpkJk  [2.9] 

k.O 

The  energy  E  is  not  an  average  because  the  correlations  are  not 
expectation  values  but  rather  geometric  defined  quantities.  The  basis 
of  the  expansion,  defined  here  for  clarity,  consists  of  the  basis 
structures  p=l  ...  m  and  the  basis  clusters  k=l  ...  M.  Therefore,  for 
a  system  of  M  equations,  there  exists  an  inverse: 

Jk  =  IHkpEp  ,  [2.10] 

p 
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provided  the  matrix  of  coefficients  5pk  is  not  singular.  Thus,  from  a 
finite  set  of  values  of  configurational  functions  a  set  of  cluster 
interactions  may  be  calculated. 

Using  Equation  2.9.  Connolly  and  Williams  (23)  calculated 
cluster  interaction  energies  from  the  first-principles  total  energies  of  a 
set  of  ordered  compounds.  In  the  present  study,  such  an  expansion 
of  the  cohesive  energies  is  used  to  examine  cluster  expansions  for 
several  lattice  types.  To  this  end.  the  excess  AF(a).  of  any 
configurational  function  F(o)  is  defined  according  to: 

AF(o)  s  F(o)  -  (1-x)  ■  F(A)  -  x-  F(B)  .  [2. 1 1) 

where  x  is  the  compositional  fraction  of  B  for  configuration  o,  and  F(A) 
and  F(B)  denote  the  properties  of  pure  A  and  B  crystals  respectively. 
Note  that  the  properties  may  have  dependence  on  variables  beside  the 
configuration  a.  For  example,  the  cohesive  energy  has  explicit 
dependence  on  the  average  atomic  volume.  By  establishing  a 
definition  of  the  excess  internal  energy,  a  reference  energy  level  for  the 
alloy  free  energy  is  defined  implicitly. 


2.3  Relaxation 

As  mentioned  above,  the  energy  depends  explicitly  on  the 
volume.  Therefore,  Equations  2.8  and  2.9  must  be  generalized  to 
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Include  explicit  dependence  of  the  cluster  interactions  of  the  volume 
V: 


Jk(V)=  IHkpEp(V)  .  (2.12) 

p 

The  correlations  themselves  depend  only  upon  the  configuration. 
However,  there  is  no  basis  to  suppose  that  the  "true"  chemical 
interactions  can  be  obtained  from  Equation  2. 1 1  using  a  single  value 
of  V.  That  is,  the  cluster  interactions  of  the  equilibrium  configuration 
of  the  disordered  alloy  probably  depend  upon  the  configurational 
energies  of  the  basis  structures  in  a  more  complex  sense  than  the 
above  relation.  If  the  equilibrium  configuration  can  be  determined  by 
constraining  each  atom  to  some  average  volume  and  then  minimizing 
the  energy  with  respect  to  this  global  volume,  then  Equation  2.11  is 
appropriate.  This  approach  is  called  global  relaxation  in  this  study. 

If  it  is  assumed  that  each  type  of  atom  "prefers"  some 
particular  volume  irrespective  of  its  environment  (location  in  a 
cluster),  then  the  interactions  can  be  relaxed  locally  —  at  the  expense 
of  excluding  non-local  effects.  This  total  relaxation  is  achieved  by 
rewriting  Equation  2. 1 1  as: 

Jk=  £  2pk  Ep(Vp) 


[2.13] 
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where  Vp  is,  for  example,  the  equilibrium  volume  of  configuration  p  at 
absolute  zero  temperature.  In  section  2.6  it  will  be  shown  that 
generally  cluster  expansions  converge  more  rapidly  with  global 
relaxation  than  with  total  relaxation. 

Ferreira  et  al  (6,  10)  "split”  the  cohesive  energy  into  chemical 
(e)  and  elastic  parts  (G).  By  assuming  that  the  equilibrium  volume  of 
an  alloy  of  fixed  concentration  is  independent  of  the  chemical  order 
they  proposed  the  so-called  e-G  representation: 

E(a,  V)  =  e(a)  +  G(x)  .  [2.141 

The  elastic  energy  G  has  explicit  dependence  on  the  concentration  x 
through  a  functional  in  which  the  equilibrium  volume  and  bulk 
modulus  of  an  alloy  are  considered  independent  of  the  state  of  order 
o.  The  chemical  energy  e  is  independent  of  volume.  This  approach  led 
to  satisfactoiy  cluster  expansions  for  semiconductor  alloys  [4),  but  the 
underlying  assumptions  —  particularly,  that  the  elastic  constants  are 
independent  of  the  state  of  order  —  are  questionable  for  intermetallics. 

Sluiter  et  al  [13J  proposed  a  simpler  scheme  that  does  not 
require  the  elastic  energy  be  independent  of  the  chemical  order.  The 
approach  consists  of  three  steps. 

1)  All  the  atoms  (A  and  B)  in  the  alloy  are  compressed  (or 
expanded)  to  the  alloy  volume  V.  Associated  with  this  compression  is 
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an  elastic  energy  called  In  this  study  the  volume  deformation  energy 
E™: 


E 


iTO  =  (l~x)(  VA  -  V)  Ba  +  x(VB-V’)  Bb 


[2.15] 


where  Va  ,  Ba  ,  Vb  .  and  Bb  are  the  equilibrium  volumes  and  bulk 
moduli  of  elements  A  and  B.  Note,  E™  is  positive. 

2)  The  atoms  are  mixed  isochorically  at  volume  Y  to  the 
equilibrium  state  of  order. 

3)  A  relaxation  energy  is  approximated: 


,  relax 


=  9^, 


12.16] 


where  W  is  the  difference  between  the  alloy  volume  V  and  the  ideal 
alloy  volume: 


W  =  V  -  (1-c)  VA  -  cVB  .  [2.17] 

and  Q  is  concentration-dependent  effective  elastic  constant.  This 
relaxation  energy  is  then  added  to  the  configurational  energy.  The 
global  relaxation  method  used  in  this  study  is  equivalent  to  steps  1 
and  2  outlined  above. 

Both  the  techniques  of  Sluiter  et  al  and  Ferreira  et  al 
described  above  are  restricted  to  volume  relaxation  (as  is  the  method 
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used  in  this  study),  i.e.,  the  atoms  are  confined  to  the  ideal  positions 
on  lattice  points.  This  treatment  ignores  the  often  significant  energy 
associated  with  relaxation  of  atoms  from  their  ideal  positions.  Zunger 
et  a l  [24]  include  two  structural  relaxation  modes  in  addition  to  the 
volume  deformation  described  above:  cell-internal  and  cell-external 
relaxations.  The  cell-external  relaxation  lowers  the  energy  by  relaxing 
the  unit-cell  vectors,  e.g.,  the  tetragonal  distortion  of  CuAu  or 
deviation  of  the  c/a  ration  in  hexagonal  structures.  Cell-internal 
relaxations  preserve  the  unit-cell  vector  lengths  but  allow  relaxation  of 
atomic  positions  not  constrained  by  the  crystal  symmetry,  e.g., 
interplanar  distances  in  CuPt-type  structures.  Note,  this  relaxation  is 
non-existent  in  certain  structures  such  as  Ll0,  Ll2.  and  D022.  These 
two  relaxation  modes  often  capture  dominant  effects  which  elude 
analyses  that  include  only  hydrostatic  relaxations.  For  example,  in 
this  study  cell-external  relaxations  lead  to  a  change  in  lattice  type  of 
the  equilibrium  structure  of  NbRu3  and  of  NbRu.  Lu  et  al  [24] 
reported  large  shifts  in  the  computed  electronic  density  of  states 
caused  by  relaxation  which  produce  better  agreement  with 
measurements.  Due  to  the  inapplicability  of  the  LMTO-ASA  to  low- 
symmetry  crystal  structures,  cell-internal  and  cell-external 
relaxations  were  not  used  extensively  in  this  study. 

Finally,  Zunger  et  al  [24]  show  that  the  effects  of  relaxation  — 
in  particular,  cell-internal  relaxation  —  are  critical  in  superlattices 
with  3  or  more  repeat  periods.  Further,  the  standard  form  of  the 
cluster  expansion,  Equation  2.8,  fails  completely  to  determine  the  cell- 
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internal  relaxation  that  produces  the  strain  energy  in  the  long-period 
superstructures. 

To  overcome  this  deficiency  Zunger,  et  al  [24]  used  a 
reciprocal- space  formulation  for  the  cluster  expansion.  The  real- 
space  interactions  J(  are  replaced  with  a  single  reciprocal-space 
function  J(k).  This  expansion  leads  to  a  rapidly  convergent  cluster 
expansion  which  predicts  accurately  the  results  of  valence  force  field 
model  of  long-period  superlattice. 


2.4  Configuration  Space 

The  method  of  extracting  cluster  interaction  energies  from 
the  configurational  energies  of  a  set  of  ordered  structures  has  been 
presented.  As  mentioned  in  the  preceding  section  a  trial  maximum 
cluster  size  is  chosen,  and  the  convergence  of  the  series  in  Equation 
2.8  must  be  established.  The  studies  cited  previously  (1.1-13)  used 
cluster  expansions  including  interactions  up  to  fifth  nearest 
neighbors.  Figures  2.1a  and  b  show  clusters  used  in  a  typical 
expansion  for  hep-  and  bcc-based  structures.  The  number  of  clusters 
chosen  determines  the  order  of  the  system  of  Equation  2.8. 

The  characteristic  functions  of  Equation  2.3  form  a  complete 
orthogonal  basis  of  configuration  space  (discrete  N-space)  [21].  The 
correlation  functions  form  a  dual  orthogonal  basis. 
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6  -  "empty"  triangle 

9  -  tetrahedron 

Figure  2. 1  b.  Typical  clusters  and  their  designated 


numbers  for  the  hep  lattice. 


The  probability  Xn(o)  of  finding  a  type-n  cluster  with 
configuration  sn  =  {i.  j . k}  (  e.g.,  {+1.  -1,  -1 . +1})  is: 

[2.18] 

Nn 

_  No.  of  n-type  clusters  with  configuration  q„ 

Total  no.  of  n-type  clusters 

The  occupation  operator  rt(p)  is  defined  as  [22]: 

r,(p)  =  [1  +  isp]  [2.19] 

=  1  if  site  p  is  occupied  by  i-type  atom,  and 
=  0  if  not. 

Using  the  occupation  operator,  the  n-type  cluster  probability  may  be 
rewritten  as  [22]: 


xjti)  =  J- £  r1(p,)rJ(pj...rk(p„)  . 

™  n  M 


[2.20] 


summing  over  all  n-type  clusters.  Combining  Equations  2.5  and  2. 19 
and  expanding  leads  to  an  expression  relating  the  cluster  correlation 
functions  with  the  cluster  probabilities  [22]: 


where  the  components  Anp  are  sums  of  products  of  the  cluster  spins 
The  requirement  that  all  cluster  probabilities  have  positive 
values,  0  <xp  <  1,  leads  to  a  set  of  constraints  on  the  correlations,  ^ 
[60]: 

1+  I  .  [2.22] 

B  >0 

not  summing  on  the  empty  cluster.  Thus,  for  any  choice  of  the 
maximum  cluster  size,  the  configuration  space  (sub-space)  is  limited 
to  a  convex  polyhedron.  The  correlation  vectors  which  form  the  K 
vertices  of  the  configurational  polyhedron,  ,  form  the  orthogonal 
basis  of  the  sub-space.  Using  barycentric  coordinates,  any  possible 
configuration,  may  be  expressed  as: 

4=  I  Pmf"  12.23] 

m  »  1 

where  the  pm  are  non-negative  coordinates.  Combining  Equation  2.23 
with  Equation  2.7  leads  to  the  configurational  description  of  the 
energy  using  clusters  p=l...K  and  structures  m=l...K: 

12241 


where 
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E<k> = I  jp4: 


tt> 

P 


p-1 


[2.25] 


is  the  energy  of  the  state  corresponding  to  vertex  k.  With  this 
representation,  Sanchez  and  de  Fontaine  [60]  have  shown  that  all 
possible  ground-state  configurations  must  correspond  to  a  vertex  or  to 
a  linear  combination  of  at  most  two  configurations  each 
corresponding  to  a  vertex  [60].  The  lattice  configurations 
corresponding  to  the  correlation  vectors  then  seem  the  natural  choice 
for  the  basis  set  of  the  expansion  of  Equation  2.8.  Indeed,  this  choice 
will  be  shown  to  be  optimum  (though  sometimes  elusive)  in  the  next 
section. 

The  enumeration  of  the  vertices  of  the  configurational 
polyhedron  is  a  difficult  problem  in  linear  programming  and  has  been 
accomplishe  only  a  limited  number  of  cluster  expansions.  The 
range  of  clus  interaction  energies  for  which  different  configurations 
are  the  ground  states  have  been  delineated  for  hep  structures  with  up 
to  second  nearest-neighbor  interactions  [61],  and  fee  and  bcc 
structures  including  up  to  fourth  nearest  neighbors  [62,  63].  These 
studies,  however,  did  not  include  many-body  interactions  and  relied 
on  the  assumption  that  nearest-  and  next  nearest-neighbor 
interactions  were  dominant.  Sanchez  and  de  Fontaine  [22]  performed 
a  full  ground-state  search  for  the  fee  lattice  including  10  sub-clusters 
of  the  joint  cluster  formed  by  the  regular  tetrahedron  (T),  composed  of 
nearest  neighbors,  and  the  regular  octahedron  (O)  which  includes  first 
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and  second  nearest  neighbors.  The  results  of  the  ground -state  search 
using  theT-O  approximation  yielded  17  possible  ordered 
configurations  whereas  the  first  and  second  nearest  neighbor 
approximation  used  by  Kanamori  [62]  yielded  only  10  of  those 
structures.  Wei  et  al  [24]  showed  that  fourth-neighbor  interactions 
are  comparable  in  magnitude  to  nearest-neighbor  interactions  for  fcc- 
based  semiconductors.  They  performed  a  ground-state  search  of  fee 
structures  including  the  many-body  interactions  while  limiting  the 
maximum  unit  cell  size  to  fifteen  atoms.  This  search,  in  contrast  to 
the  vertex  enumeration  studies,  calculates  the  correlation  vectors  of 
all  possible  independent  structures  and  then,  using  the  cluster 
expansion  equation,  the  cohesive  energies.  Such  a  method  is 
employed  in  this  study,  in  a  limited  manner,  and  is  discussed  in 
Section  2.6. 

2.8  How  to  choose  a  basis  set 

For  the  alloy  systems  in  this  study  (as  well  as  most  other 
studies)  interaction  ranges  longer  than  nearest  neighbors  are 
required.  Without  knowing  the  vertices  of  the  configurational 
polyhedron,  a  means  of  estimating  the  optimum  choice  of  basis 
structures  for  a  given  maximum  cluster  size  is  required.  Three 
techniques  are  considered  --  yielding  different  conclusions: 

(i)  Random  alloy  convergence. 
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(ill  Minimum  eigenvalue  norm. 

(ill)  Maximum  determinant  method. 

The  expectation  value  of  the  point  cluster,  E,i.  depends  only 
upon  the  alloy  concentration  c: 

=  2c  -  1  .  [2.261 

In  a  random  alloy  each  lattice  site  is  equivalent,  thus  each  site  in  a 
cluster  has  the  same  value  of  £i ,  and  the  correlation  function  of  a 
cluster  with  k  sites 

random  l 

Sr  =  (2c  -  1)  .  [2.27] 

Notably,  at  c  =  0.5  .  all  terms  vanish  except  at  the  k  =  0  term,  i.e.,  the 
empty  cluster: 

— random 

F  (c  =  0.5}  =  J0  .  [2.28] 

— random 

If  the  value  of  F  can  be  calculated  directly,  using  such  techniques 
as  the  CPA-GPM,  then  Equation  2.21  can  be  used  to  evaluate  the 
accuracy  of  the  expansion. 

Combining  Equations  2.8  and  2.9: 

°F  =  I  %  Hk'p  Fp  . 

k.  p 


[2.29] 
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=  I  °QpFP 

P 

where, 

°QP  =  I  nkHk'p.  [2.30] 

k 

Thus,  the  duster  expanded  configurational  property  of  a  random  alloy 
is: 


rand. 

<U 


rand. 


[2.31] 


where. 


rand.  _  m  , „  k  ** 

Qp  =I(2c-l)  Hkp  .  [2.32] 

k 

The  sensitivity  of  the  random  alloy  property  is: 

— rand. 

12331 

Ferreira  et  al  [12]  use  this  result  to  determine  the  stability  of  the 
expansion  of  the  excess  cohesive  energy.  They  integrate  the  norm  of 

rind. 

Qp  over  all  concentrations: 


„  .  f  * _ rand.  _nnd.  , 

y(k.p)=IjoQp  QP  dc 


[2.34] 
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giving: 


y(k,p)  =  II  Hkp  Hp,  x 

k.  I  P 


1  -  (-1  ) 


k  +  U  1 


2(k  +  I  +  1) 


[2.35j 


The  set  of  basis  structures.  k=l,...n,  that  yields  the  minimum  value  of 
y  corresponds  to  the  expansion  least  sensitive  to  errors  in  the 
cohesive  energy  for  the  choice  of  clusters,  p=l,  ...  n. 

The  second  and  third  methods  also  use  estimates  of  the 
sensitivity  of  the  interactions  to  errors  in  the  excess  cohesive  energies. 
From  Equation  2.9,  the  sensitivity  of  Jk  is: 

5J  «  S  *  8F  .  [2.361 


Using  a  standard  result  of  numerical  analysis  (see  Wilkinson  [64]). 
Equation  2.35  may  be  refined  to: 

=  CIS1)  Igt  .  [2.37] 

w 

where  C(A)  is  the  condition  number  of  A,  defined  by: 


M 

M 


C(A)=  (A||A'| 


[2.38] 
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and  |A|  denotes  some  type  of  matrix  norm.  The  norm  used  in  this 
study  is: 

T  1/2 

|A|  =  (max.  eigenvalue  of  A  A)  [2.39] 

-i 

The  choice  of  E  which  minimizes  C(E  )  is  deemed  the  optimum  basis 
structure. 

Finally,  noting  that: 

5Jk  =  I  Ekp  8F  .  (2.401 

p 

the  uncertainty  in  each  potential  5Jk  depends  upon  the  sum  of  the 

-> 

absolute  values  of  each  element  in  the  corresponding  row  of  Ekp  .  By 
Kramer's  rule: 

Det  (E,j)| 

-kp  =  r  i  *  k.  P  •  [2.41] 

Det  (S)| 

it  is  clear  that  the  condition  number  of  Equation  2.36  is  minimized  for 
the  choice  of  S  for  which  the  determinant  is  a  maximum  and  the 
determinant  of  each  minor  is  minimum.  Thus,  the  simple  norm,  the 
determinant  of  the  correlation  matrix,  E  is  adopted. 

Figures  1.6  a  and  b  show  the  bcc-  and  hep-based  structures 
considered  as  candidates  for  use  as  basis  structures.  They  are 
chosen  due  to  their  small  unit-cell  sizes,  a  necessity  for  the  LDA 


Table  2.1  a.  Correlation  functions  for  selected  bcc  structures  of 
Figure  1.1  a. 
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Structure 

No. 

O 

• 

/ 

«=* 

u 

♦ 

bcc  (A) 

1 

1 

1 

1 

1 

1 

1 

D03  (AgB) 

2 

1 

1/2 

0 

0 

-1/2 

-1 

B2(AB) 

3 

1 

0 

-1 

1 

0 

1 

B32  1AB) 

H 

1 

0 

0 

-1 

0 

1 

5 

1 

-1/2 

0 

0 

1/2 

-1 

bcc  (B) 

6 

1 

-1 

1 

1 

-1 

1 

T1  (AB) 

n 

1 

0 

0 

-1/3 

0 

-1/3 

T2  (AB) 

8 

1 

0 

-1/4 

1/3 

0 

0 

T3  (AB) 

9 

1 

0 

0 

0 

0 

-1/3 

T4  (AB) 

10 

1 

0 

0 

1/3 

0 

-1/3 

T5  (A3B) 

11 

1 

1/2 

0 

2/3 

1/6 

1/3 

T5(ABa) 

12 

1 

-1/2 

0 

2/3 

-1/6 

1/3 

1 

1 

4 

3 

12 

6 
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Table  2.1  b.  Correlation  functions  for  selected  hep  structures  of 

Figure  1.1  b. 


Structure 


hep  (A)  1 


DOi9  (A3B)  2 


Pmma  (AB)  3 


Pmmn  (AB) 


DO, 9  (AB3)  5 


hep  (B)  6 


LlRh(AB) 


HI  (AB)  8 


Degenera 


1 

1 

1 

1 

1 

1 

1 

1 

1/2 

0 

0 

1 

-1/2 

-1/2 

1 

0 

-1/3 

-1/3 

1 

0 

0 

1 

0 

1/3 

-1/3 

-1 

0 

0 

1 

-1/2 

0 

0 

1 

1/2 

1/2 

1 

-1 

1 

1 

1 

-1 

-1 

1 

0 

-1 

1 

-1 

0 

0 

1 

0 

0 

-1/3 

0 

0 

0 

1 

1 

3 

3 

3 

6 

1 
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Table  2. 1  c.  Correlation  functions  for  selected 
hep  structures.  Isotropic  interactions  are 
assumed. 


Structure 

No. 

0 

• 

/ 

•=» 

Q 

H 

hep  (A) 

1 

1 

1 

1 

1 

l 

1 

DO  is  (A3B) 

2 

1 

1/2 

0 

1 

-1/2 

-1 

Pmma  (AB) 

3 

1 

0 

-1/3 

1 

0 

1 

Pmmn  (AB) 

H 

1 

0 

0 

-1 

0 

-1 

5 

1 

-1/2 

0 

1 

1/2 

-1 

hep  (B) 

6 

1 

-1 

1 

-1 

-1 

1 

LIRh  (AB) 

n 

1 

0 

0 

-1 

0 

-1 

HI  (AB) 

8 

1 

0 

-1/6 

0 

0 

0 

n 

1 

1 

6 

3 

8 

2 
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calculations.  The  correlation  functions  of  clusters  of  up  to  fifth 
nearest  neighbors  are  given  for  each  of  these  structures  in  Tables  2. 1 
a.  b,  and  c.  The  first  column  is  the  crystal  structure,  and  the  second 
column  is  the  structure  number  provided  for  convenience.  The  rest  of 
the  columns  are  the  correlation  functions  for  each  of  the  cluster  types 
represented  by  the  column  heading.  The  final  row  is  the  cluster 
degeneracy. 

Table  2.2  lists  the  best  choices  for  a  basis  set  according  to 
criteria  (i)  maximum  determinant,  (ii)  smallest  eigenvalue,  and 
(iii)  random  alloy  parameter  for  each  of  three  different  sets  of  basis 
clusters  A.  B,  and  C.  For  cluster  set  A,  the  basis-structure  set  1  -6  is 
the  optimum  choice  for  each  of  the  criteria.  This  is  not  surprising 
since  the  correlation  vectors  of  these  structures  are  the  vertices  of  the 
associated  configuration  polyhedron.  For  cluster  set  B.  the  three 
basis  sets  (1-61  +  (7,  9,  or  10)  yield  the  same  determinant.  This 
degeneracy  is  removed  by  criteria  (ii)  and  (iii).  Notably,  criterion  (iii) 
predicts  an  optimum  basis  set  which  excludes  structure  3  (CsCl 
structure)  for  both  cluster  bases  B  and  C  whereas  criterion  (i)  is 
optimized  by  basis  sets  which  include  the  structure.  These 
convergence  criteria  for  the  expansion  basis  are  used  in  the  next 
section  for  numerical  evaluation  of  the  convergence  of  the  cluster 
expansion  for  the  three  different  systems. 
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Table  2.2.  Optimum  choices  of  basis  sets  for  bcc  structures  for  each  of  the  three 
criteria  described  in  section  2.5.  Note:  the  structure  numbers  in  the  basis  sets  are 
given  in  Table  2.1  a. 


Basis  clusters* 

Optimization  criterion 

number  of 

determinant 

eigenvalues 

random  alloy 

sub-clusters 

A 

[1-6] 

[1-6] 

[1-6] 

6 

B 

[l-6]  +  9 

[1-61  +  9 

1245678 

7 

[l-6]  +  7 

[1-61  +  7 

1245689 

[1-61+10 

C 

[1-7J  +  10 

1  245  689  10 

1  245678  10 

8 

[1-71  +  8 

12456789 

1  24  5  67  1012 

Cluster  set  A  consists  of  clusters  *  1, 2, 3, 4, 5,  and  6  [see  Figure  2.1a]. 


"  B . 1,2, 3, 4, 5, 6,  and  7. 


C  "  “ 


1, 2, 3, 4, 5, 6, 7,  and  9. 


Recall  that  Equation  2.8  allows  an  estimate  of  any 
configurational  property  for  a  given  choice  of  a  basis  set  of  structures 
and  maximum  cluster  size.  The  obvious  method  of  gauging 
convergence  of  the  cluster  expansion  is  to  evaluate  the  expansion 
coefficients,  l.e.,  the  cluster  interactions.  The  interactions  should 
vanish  with  interaction  distance  and  with  cluster  complexity.  The 
value  of  each  type  of  interaction  should  also  converge  to  a  stable  value 
as  the  size  of  the  basis  set  grows.  Thus,  any  configurational  property, 
for  example,  the  cohesive  energy  calculated  using  Equation  2.8, 
should  converge  to  a  value  in  a  stable  manner;  note  though,  this  is 
not  necessarily  uniform  convergence.  This  second  observation  leads 
to  a  direct  method  for  evaluating  the  accuracy  of  the  cluster 
expansion:  comparing  the  value  of  a  configurational  property 
evaluated  directly  with  the  cluster-expanded  value.  In  this  section, 
the  convergence  of  a  cluster  expansion  is  judged  by  estimating  the 
cohesive  energy  of  an  ordered  compound  using  Equation  2.8  and 
comparing  the  value  with  the  LDA  value  from  Table  1.2.  These 
assessments  are  performed  for  bcc-  and  hep-based  lattice  structures 
in  Zr-Nb,  Nb-Ru,  and  Zr-Ru.  The  success  with  which  the  cohesive 
energies  of  the  ordered  compounds  are  predicted  lends  confidence  to 
the  predictions  for  disordei  ed  alloys. 
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Zr-Nb 

Figure  2.2  shows  the  ground-state  excess  cohesive  energies 
of  selected  hep-  and  bcc-based  compounds  in  the  Zr-Nb  system. 
Clearly  there  are  no  stable  intermediate  compounds  --  as  expected  in 
system  with  a  predominant  miscibility  gap.  Using  the  different  basis 
sets  advocated  by  each  of  the  three  optimization  techniques  shown  in 
Table  2.2,  the  excess  cohesive  energies  of  extra  structures  were 
calculated  using  the  cluster  expansion  with  both  global  and  total 
volume  relaxation.  The  resultant  errors  for  each  expansion  are 
summarized  in  Table  2.3  using  different  error  norms.  The  expansion 
using  clusters  up  to  second  nearest  neighbors  converges  quickly  for 
the  globally  relaxed  case,  but  at  least  third  nearest  neighbors  are 
needed  for  the  totally  relaxed  expansion.  Omitting  the  B2  structure 
from  the  basis  set  for  expansions  of  up  to  fifth  nearest  neighbors 
yields  predictions  with  a  larger  error  norm  than  the  second  neighbor 
expansion.  This  effect  is  due,  in  part,  to  the  strong  dependence  of  the 
B2  energy  on  the  nearest  neighbor  interaction.  Figure  2.3  shows  the 
unstable  description  of  this  interaction  for  expansions  that  omit  the 
B2  structure  from  the  basis(curves  #3,  4.  and  5)  and  the  very  stable 
result  of  basis  sets  of  varying  size  that  include  B2  (curves  #1-3). 

Figure  2.4  shows  the  derived  interactions  for  the  "best" 
expansions  with  seven  (up  to  third  nearest  neighbors)  and  eight 
cluster  interactions  (up  to  fifth  nearest  neighbors).  For  each  of  the 
expansions  the  nearest  neighbor  interactions  are  nearly  identical. 


Table  2.3.  Errors  of  cluster  expanded  excess  structural  energy  of  "extra"  bcc 
structures  (i.e.,  those  not  included  in  basis  set)  in  Zr-Nb.  Values  are  given  for 
globally  relaxed  potentials  and  (totally  relaxed). 
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cluster 

basis. 

number 

basis 

max  lerrorl 

ImRy.l 

£  II  (error)2] 1/2 

[mRy.J 

1  ,  Jerroif  1  1/2 

[percent] 

A  6 

1-6 

.59  (2.52) 

.19  (.91) 

4.0(21.0) 

B  7 

1-7 

.49(1.23) 

.20  (.35) 

4.1  (7.4) 

B  7 

1-6  +  9 

.32  (.69) 

.13  (.22) 

3.0  (5.4) 

B  7 

1245678 

.89(1.22) 

.32(41) 

5.5  (7.7) 

B  7 

1245689 

.64  (.69) 

.23  (.22) 

4.1  (54) 

C8 

1-8 

.13(42) 

.08  (.11) 

1.9  (2.3) 

C  8 

12456789 

.53  (.52) 

.27  (.14) 

4.6  (2.3) 

1.0 


Nb  Zr 

Concentration  of  Zr 

Figure  2.3  Nearest  neighbor  interactions  calculated  using  cluster 
expansions  with  different  bases  and  global  relaxation. 
1*123456 

2- 1234569 

3- 12345678 

4- 12456789 

5- 1245689 

6- 1245678 
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The  second  nearest-neighbor  interactions  are  identical  for  the 
expansions  including  second  nearest  and  third  nearest  neighbors,  but 
the  inclusion  of  fifth  nearest  neighbors  causes  a  shift  in  the 
interaction.  The  third-  and  fifth-neighbor  interactions  are  small 
(approximately  0.05  mRy./atom)  and  nearly  constant  over  the  entire 
concentration  range.  The  magnitudes  of  the  pair  interactions  are 
closely  related  to  the  relative  stabilities  of  alloys  of  equi-atomic 
concentrations.  At  a  concentration  of  50%,  the  second  nearest- 
neighbor  interaction  is  approximately  2.5  times  as  large  as  the 
nearest-neighbor  interaction.  The  correlation  functions  of  the  two 
pair  clusters  for  the  B2  and  B32  structures  are  =-1,^=1, 

5“2=  5**  =  -1  (see  Table  2. 1).  Thus,  the  contribution  to  the  B2 

formation  energy  of  the  second  nearest  neighbor  interaction  is  positive 
(ferromagnetic)  and  destabilizing.  However,  the  B32  structure  is 
stabilized  by  the  relatively  large  negative  (-2.8  mRy/atom)  second 
neighbor-interaction,  and  the  B32  compound  is  much  more  stable 
than  the  B2.  Why  is  the  the  B32  structure  not  stable?  The  answer 
lies  in  the  empty  cluster  interaction  which  can  be  associated  with  the 
random  alloy  formation  energy  .  The  random  formation  energies  (both 
globally  and  totally  relaxed)  as  a  function  of  concentration  with  those 
of  the  bcc-based  compounds  are  shown  in  Figure  2.5.  Recall  from 
Equation  2.26  that  the  formation  energy  for  a  random  alloy  with  50% 
concentration  is  equal  to  the  value  of  the  effective  interaction  Jo.  for 
the  empty  cluster.  As  seen  in  Figure  2.5  the  value  of  Jo  is 
approximately  4. 1  mRy/atom.  This  interaction  is  the  dominant  and 
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destabilizing  contribution  to  the  formation  energies  of  the  compounds 
near  equi-atomic  concentrations.  The  metastable,  intermediate  bcc- 
based  compounds  2,  4.  5.  7,  and  9  are  stable  with  respect  to  the 
random  alloy,  which  suggests  short-range  ordering.  The  difference 
between  the  formation  energy  of  a  compound  and  the  random 
formation  energy  at  that  concentration  is  the  ordering  energy.  Thus, 
the  unstable  B32  compound  has  a  negative  ordering  energy  of  2.5 
mRy/atom. 

A  metastable  hexagonal  to  phase  varying  from  P6/mmm 
space  group  to  P3ml  with  increasing  Nb  content  may  be  retained  by 
quenching  165).  The  hexagonal  phases  1-6  are  unstable  with  respect 
to  the  random  bcc  alloy,  however,  so  no  connection  is  made  between 
the  results  of  this  study  and  those  of  [651. 

The  cluster  expansion  on  the  hexagonal  lattice  requires 
tracing  the  effects  of  clusters  that  lie  completely  within  the  basal 
plane  or  have  out-of-plane  vertices.  As  seen  in  Table  2. 1  b,  the 
correlation  functions  depend  upon  the  orientation  of  otherwise 
identical  clusters.  If  isotropic  interactions  are  assumed,  then  the  first 
neighbor  and  first-neighbor  triangles  can  be  combined,  yielding 
effective  correlation  functions  (see  Table  2.1  c).  Such  an  assumption 
renders  the  correlation  vectors  for  structures  4  (Pmmn)  and  7  (LiRh) 
identical  up  to  (and  including)  second  nearest  neighbors.  This 
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Concentration  of  Zr 


Figure  2.4.  Pair  interactions  in  Zr-Nb  derived  from  cluster  expansions 
of  different  range.  The  size  of  the  basis  set  is  indicated  by  the  symbol 
type: 

a 

■  Nearest  neighbors  only.  Basis  =  1  2  3  4  5  6 
O  3rd  nearest  neighbors  included.  Basis  =  1  2  3  4  5  6  7 

- 5  nearest  neighbors  included.  Basis  =12345678 


Excess  structural  energy  [mRy/atom] 
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Nb  Zr 


Concentration  of  Zr 

Figure  2.5.  Random  energy  for  bcc  alloys  in  Zr-Nb  calculated  using 
both  globally  and  totally  relaxed  cluster  expansions.  Superimposed 
are  the  excess  structural  energies  of  bcc  compounds  from  Figure  2.4 
using  Zr-bcc  and  Nb-bcc  as  reference  energies. 
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observation  suggests  an  obvious  check  for  the  isotropy  assumption:  if 
the  formation  energies  for  compounds  with  crystal  structures  4  and  7 
are  approximately  equal,  then  the  assumption  of  isotropic  short-range 
interactions  is  valid.  This  is  the  case  for  the  Zr-Nb  system.  Table  1.6 
b  and  Figure  2.6  show  the  cohesive  energies  of  the  two  compounds 
(curves  #2  and  #3}  as  a  function  of  average  atomic  volume  which  are 
nearly  identical  (AE=0.1  mRy./atom).  Curve  #1  is  the  projected 
cohesive  energy  for  both  compounds  using  a  cluster  expansion  of  up 
to  nearest  neighbors  and  assuming  isotropy.  Clearly,  the  error,  1.5 
mRy./atom.  indicates  an  unconverged  expansion.  If  second  nearest 
neighbors  are  included  in  the  expansion,  and  structure  4  is  included 
in  the  basis,  then  curve  #2  is  the  projected  energy  for  structure  #7. 
Thus,  for  Zr-Nb  the  hexagonal  cluster  interactions  are  reasonably 
converged  using  the  isotropy  assumption  and  including  clusters  of  up 
to  second  nearest  neighbors. 

The  excess  volumes  and  bulk  moduli  of  the  hep-,  and  bcc- 
based  compounds  were  cluster  expanded  using  the  total  relaxation. 
The  effective  interactions  are  shown  in  Table  2.4  for  each  type  of 
cluster  used  in  the  expansion.  The  bcc  expansion  includes  2nd 
nearest  neighbors.  The  hep  expansion  includes  nearest-neighbor  pair 
interactions  and  includes  the  full  isotropy  assumptions  made  in  the 
energy  expansion.  The  convergence  of  the  bulk  modulus  effective 
interactions  Bj  is  rapid.  The  four-body  interaction  is  four  to  five 
orders  of  magnitude  below  1  Mbar,  the  nominal  size  of  the  bulk 
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Average  Wigner-Seitz  radius  [a.u.] 


Figure  2.6.  Excess  structural  energies  of  ZrNb  with  Pmmn  and  LiRh 
structures  from  fitted  Morse  potentials  and  cluster  expansion. 

The  cluster  expansion  includes  only  nearest  neighbors  and  uses 
isotropy  assumption  and  is  indicated  by  curve  #1 .  Curve  #3  is  the 
Morse  fit  for  the  LiRh  structure.  Curve  #2  is  the  the  fit  for  the  Pmmn 
structure.  When  second  nearest  neighbors  are  included  in  the 
expansion,  and  Pmmn  is  included  in  the  basis  set,  curve  #2  becomes 
the  projected  energy  for  the  LiRh  structure  since  the  correlation  vectors 
for  LiRh  and  Pmmn  are  identical  when  isotropic  interactions  are 
assumed. 
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moduli  in  this  sytem.  The  convergence  of  the  volume  expansion  is 
rapid  for  the  bcc-based  structures,  and  sufficient  in  the  hep  alloys. 

Table  2.4.  The  interaction  coefficients  for  cluster  expansion  of  bulk  modulus  and 
the  volume  of  (a)  bcc  structures  and  (b)  hep  structures  in  Zr-Nb.  The  cluster  labels 
are  from  Figure  2. 1 . 

(a)  _  _ 


Cluster 

Effective  interactions 

(label) 

Bj  [Mbar] 

Vi  [a.u.3] 

Empty  (1) 

-0.049120 

-0.937820 

Point  (2) 

0.001500 

-0.069780 

Pair  (3) 

0.009875 

0.114747 

2n<*  pair  (4) 

0.003083 

0.144064 

Triangle  (5) 

-0.000120 

0.005817 

Tetrahedron  (6) 

0.000062 

0.007776 

Cluster 

Effective  interactions 

(label) 

Bi  [Mbar) 

Vj  [a.u.3) 

Empty  (1) 

-0.058440 

-0.297440 

Point  (2) 

-0.000250 

0.024432 

Pair  (3  4) 

0.009563 

0.061917 

Triangle  (6  7  8) 

0.000031 

-0.003060 

Tetrahedron  (9) 

0.000531 

-0.037020 
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Nb-Ru 

The  ground-state  energies  of  selected  compounds  in  the  Nb- 
Ru  system  are  shown  in  Figure  2.7.  The  stable-equilibrium  structure 
of  the  NbRu3  compound  is  predicted  to  be  DO  19,  but  the  margin  of 
stability  with  respect  to  the  LI  2  structure  (0.1  mRy./atom)  is  on  the 
order  of  the  precision  of  the  total  energy  calculations.  This  delicate 
balance  provokes  the  question  of  whether  a  different  fee  or  hexagonal 
structure  might  have  an  even  more  negative  (stable)  formation  energy. 

To  resolve  this  question  the  cohesive  energies  of  several 
candidate  structures  were  predicted  using  the  cluster  expansion.  For 
the  fee  and  hexagonal  cluster  expansions,  interactions  of  up  to  second 
nearest  neighbors  were  included.  The  correlation  functions  of  all  fcc- 
and  hep-based  structures  with  stoichiometry  AB3  and  up  to  16  atoms 
per  unit  cell  were  calculated.  The  unit  cells  were  restricted  to  those 
occupying  a  supercell  obtained  by  doubling  the  standard  unit  cell  size 
in  each  principal  direction. 

The  cluster  interactions  were  obtained  with  a  globally  relaxed 
expansion  using  the  basis  set  deemed  optimum  by  the  maximum 
determinant  technique.  The  assumption  of  isotropy  in  the  cluster 
expansion  of  hexagonal  lattice  structures  is  invalid  as  shown  in  Figure 
2.8  a.  since  the  cohesive  energies  of  the  compounds  with  structures  4 
and  7  have  a  significant  difference  (1.8  mlty/atom).  The  projected 


Formation  energy  [mRy/atom] 


Figure  2.7.  Ground-state  energies  of  selected  compounds  in  Nb-Ru, 
and  random  energies  of  hep  and  bcc  alloys. 

-  Bcc  random  alloy  (Totally  relaxed) 

- "  "  "  (Globally  relaxed) 

.  Hep  random  alloy  (T otally  relaxed) 

- "  "  "  (Globally  relaxed) 


Structural  energy  (mRy.]  Structural  energy  [mRy.]  Structural  energy  [mRy.] 
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Average  Wigner-Seitz  radius  [a.u.] 


Figure  2.8  a.  Projected  and  calculated 
structural  energies  for  NbRu  Pmmn- 
and  LiRh  structures.  Expansion 
includes  nearest  neighbors,  and 
isotropic  pair  interactions  are  assumed. 
Note:  projected  energy  for  LiRh  and  Pmmn 
structures  are  identical  for  this  expansion. 


Figure  2.8  b.  Projected  and  calculated 
structural  energies  for  NbRu  Pmmn,  LiRh, 
and  HI  structures.  Expansion 
includes  nearest  neighbors,  and  anisotropic 
pair  interactions  are  assumed. 


Figure  2.8  c.  Projected  and  calculated 
structural  energies  for  NbRu  with  Hi 
structure.  Expansion  includes  second 
nearest  neighbors  with  anisotropic  pair 
interactions. 
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energy  for  the  two  structures  using  a  nearest-neighbors  expansion 
and  isotropic  interactions  is  also  shown  in  Figure  2.8  a  with  an 
obviously  large  error  (~  5  mRy./atom).  Figure  2.8  b  shows  the 
projected  energies  for  structures  7  and  8  in  a  nearest-neighbors 
expansion  with  anisotropic  pair  interactions.  The  discrepancies  are 
as  large  as  15  mRy./atom.  Including  second  nearest  neighbors 
improves  the  expansion  significantly.  The  projected  energy  for 
structure  8  using  this  expansion  is  shown  in  Figure  2.8  c.  The  pair 
interactions  of  this  well  converged  expansion  are  shown  Figure  2.9. 

No  hexagonal  or  cubic  structures  were  found  with  projected 
formation  energies  more  negative  than  that  of  the  DO  19.  However,  it 
must  be  noted  that  the  LMTO-ASA  approximation  loses  validity  for 
hexagonal  structures  varying  significantly  from  ideal  close  packing. 
Thus,  an  improved  energy  analysis  will  require  full-potential  total- 
energy  calculations. 

The  Nb-rlch  DO3  is  also  a  predicted  ground  state,  however 
no  observation  of  such  a  phase  has  been  reported.  A  ground-state 
search  similar  to  the  one  for  NbRu3  was  performed  yielding  no 
compound  more  stable  than  DO3. 

Table  2.4  shows  the  results  of  a  convergence  study  of  the  bcc 
expansion.  The  last  column,  the  percent  error  norm,  is  of  comparable 
magnitude  to  the  Zr-Nb  bcc  cluster  expansion.  Figure  2.10  shows  the 
cohesive  energy  of  the  random  bcc  alloys  calculated  using  global 
relaxation  and  different  basis  sets.  The  random  energy  appears  to  be 
well  converged  using  the  basis  set  {1  2  3  4  5  6  7).  The 
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Figure  2.9.  Pair  interactions  derived  from  cluster  expansion 
of  Nb-Ru  hep  structural  energies. 


Table  2.5.  Errors  of  cluster  expanded  excess  structural  energy  of 
"extra"  bcc  structures  (l.e. .  those  not  included  in  basis  set)  in  Nb-Ru. 


cluster 

basis. 

number 


basis 


1245678 


1-6  +  9 


1245689 


1  24  56  7  8  10 


1-7+10 


1  24  56  8  9  10 


1 2456789 


1-7  +  8 


1  24  56  7  10  12 


max  1  error  1 

[mRy.l 

4-  lI|enror|a] 1/2 

*1  n 

Imlty.) 

8.5  (5.8) 

2.8  (1.5) 

4.2  (3.5) 

1.0  (1.2) 

6.6  (10.4) 

1.8  (2.4) 

4.6  (3.6) 

1.0  (.98) 

5.0  (4.6) 

1.7  (1.3) 

8.5  (8.5) 

2.6  (2.6) 

4.0  (3.6) 

1.2  (1.1) 

5.6  (5.6) 

2.1  (2.0) 

5.9  (5.7) 

1.8  (1.8) 

5.9  (5.7) 

1.6  (1.8) 

5.0  (4.7) 

1.7  (1.7) 

10.9  (5.5) 


4.1  (5.2) 


6.4  (7.9) 


4.0  (3.7) 


6.1  (4.8) 


8.9  (8.7) 


4.8  (4.6) 


7.2  (6.9) 


7.0  (6.7) 


7.5  (7.2) 


6.0  (6.1) 


Structural  energy  [mRyVatom] 


Figure  2.10.  Excess  structural  energy  of  random  bcc  alloys  in 
Nb-Ru  derived  from  globally  relaxed  cluster  expansions  of 
different  approximations.  Symbols  denote  formation  energies  of 
selected  bcc  compounds.  The  different  cluster  expansions  are: 

1  -  2nd  neighbors,  basis  =  1  2  3  4  5  6  A 

2 -  3rd neighbors,  basis  =  1245678  B 

3- "  "  "  =1 234567  B 

4  -  3rd  neighbors  (B)+  linear  triplet,  basis  =  123456912 

5 -  51h neighbors,  basis  =  123456710  C 

6- "  "  "  =  1  245678  10  C 
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cluster  Interactions  obtained  from  the  expansion  with  a  basis  set  {1  2 
3  4  5  6  7  8}  and  interactions  of  up  to  fifth  nearest  neighbors  are 
shown  in  Figure  2. 1 1. 

A  NbRu  compound  with  the  Llo  structure  is  predicted  to  be 
stable  as  expected.  Figure  2.12  shows  the  strong  dependence  of  the 
cohesive  energy  on  the  tetragonality  of  the  phase.  The  calculated  c/a 
ratio  is  1.122,  which  is  identical  to  the  measured  value  reported  by 
Hurley  and  Brophy  (49].  Das  et  al  (38]  observed  a  tetragonal  to 
orthorhombic  transition  at  46  %  Ru  (see  Figure  1.7),  and  the  average 
of  the  orthorhombic  a/c  and  b/c  ratios  measured  at  a  concentration 
of  51%  Ru  is  1. 126.  An  explanation  of  the  stabilization  of  the  NbRu 
phase  by  a  tetragonal  distortion  lies  in  an  analysis  of  the  pair 
interactions.  Figure  2.13  displays  the  interactions  derived  from  a 
globally  relaxed  cluster  expansion  at  the  optimum  atomic  volumes  for 
each  of  three  c/a  ratios  --  corresponding  to  a  B2  phase,  Llo  phase 
with  c/a=1.122,  and  a  cubic  Llo  phase.  The  Llo  phase  is  stabilized 
by  a  change  in  sign  of  the  third  nearest  neighbor  interaction  in 
conjunction  with  a  slight  increase  in  the  nearest  neighbor  interaction. 
Further  tetragonal  distortion  is  destabilized  by  the  "mixing"  of  what 
were  the  first  and  second  nearest  neighbors  of  a  bcc  lattice  to  nearest 
neighbors  in  an  fee  lattice,  i.e.,  the  cubic  Llo  has  a  correlation 
function  of  zero  for  the  nearest  neighbor  cluster. 


Cluster  interaction  (mRy./atom) 
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Figure  2.11.  Cluster  interactions  of  bcc  alloys  in  Nb-Ru  including 
5th  nearest  neighbor  interactions  and  global  relaxation. 

- 5th  nearest  neighbors 

. Tetrahedron 

- Triangle 

- 3rd  nearest  neighbors 

- 2nd  nearest  neighbors 

- Nearest  neighbors 


Structural  energy  [mRy/atom] 


Figure  2.12.  Structural  energy  of  NbRu-L10  as  a  function  of  the 
c/a  ratio.  At  c/a=1  the  structure  is  equivalent  to  B2,  and  at 
c/a  =  1 .414  the  structure  is  a  cubic  LI  0.  The  dashed  line  indicates 
the  energy  corresponding  to  the  line  of  stability  between  Nb3Ru-D03 
and  NbRu3-D019. 


lir  interactions  (mRy/atom) 


F'n,,re  213  Pair  interactions  as  a  function  of  the  c/a  ratio.  Note 
that  the  3*  nearest  neighbor  interactions  for  the  foe  lathee  are 

forced  to  zero. 


Table  2.6  shows  the  interaction  coefficients  from  a  cluster 
expansion  of  the  excess  bulk  moduli  and  volumes  of  bcc-,  hep-,  and 
fcc-based  structures  in  Nb-Ru.  The  convergence  of  the  bulk  modulus 
expansion  appears  to  be  more  rapid  than  for  the  volume  expansion. 
The  tetrahedron  interaction  for  the  hep  volume  expansion  is  relatively 
large  (one  order  of  magnitude  less  than  the  pair  interaction)  as  is  the 
case  with  Zr-Nb. 
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Table  2.6.  The  Interaction  coefficients  for  cluster  expansion  of  bulk  modulus  and 
the  volume  of  (a)  bcc  structures  (b)  hep  structures,  and  (c)  fee  structures  in  Nb-Ru. 
The  cluster  labels  are  from  Figure  2.1. 

(a)  _ 


Cluster 

Effective  interactions  | 

(label) 

Bj  [mbar] 

Vj  [a.u.3] 

Empty  (1) 

0.010750 

-1.502580 

Point  (2) 

0.027500 

-0.640950 

Pair  (3) 

-0.007250 

0.185959 

2nd  pair  (4) 

0.007500 

0.223879 

Triangle  (5) 

-0.002290 

0.053412 

Tetrahedron  (6) 

-0.000710 

0.014518 

(b) 


Cluster 

Effective  interactions  j 

(label) 

Bf  (mbar) 

Vj  [a.u.3] 

Empty  (1) 

-0.072940 

-1.262280 

Point  (2) 

-0.004750 

-0.041540 

Pair  (3  4) 

0.011813 

0.199157 

Triangle  (6  7  8) 

0.000594 

0.005192 

Tetrahedron  (9) 

0.001031 

0.033671 

Cluster 

Effective  Interactions  f 

Bi  [mbar] 

Vj  [a.u.3] 

Empty 

-0.039790 

-1.607730 

Point 

-0.018030 

0.116430 

Pair 

0.008256 

0.255239 

2nd  pair 

0.002253 

-0.014550 

Triangle 

0.001128 

0.020825 

Tetrahedron 

-0.004000 

0.011547 
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Zr-Ru 

Table  2.7  shows  the  error  norms  for  bcc  cluster  expansions  using  the 
basis  sets  optimized  according  to  criterion  i.  The  expansion  seems  to 
converge  rapidly  with  the  inclusion  of  3rd  nearest  neighbors,  but 
including  5th  nearest  neighbors  does  not  improve  the  expansion. 
Figure  2.14  shows  the  random  alloy  formation  energies  for  hep  and 
bcc  alloys  in  Zr-Ru.  The  energies  were  calculated  using  both  total  and 
global  relaxation.  The  difference  between  these  two  extreme  methods 
of  volume,  approximately  7  mRy/atom  for  both  lattice  types,  is 
significant,  indicating  the  large  effect  relaxation  has  in  this  system. 
The  ordering  energy  of  ZrRu  with  the  B2  structure  is  38  mRy/atom. 
Based  on  the  total-energy  calculations  reported  by  Lu  et  al  [24]  the 
ordering  energies  for  NiAl  (B2)  and  Ni3Al  are  50  mRy/atom  and  13.5 
mRy/atom  respectively. 


Table  2.7.  Errors  of  cluster  expanded  excess  structural  energy  of 
"extra"  bcc  structures  (i.e.,  those  not  included  in  basis  set)  in  Zr-Ru. 
Values  are  given  for  globally  relaxed  potentials  and  (totally  relaxed). 
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cluster 

basis. 

number 

basis 

max  1  error  1 

ImRy.) 

±  lX|error|2] 1/2 

ii  n 

(mRy.l 

1  ,  Jerrorf,  1/2 

[percent) 

A  6 

1-6 

12.3  (6.5) 

3.37  (1.7) 

9.7  (5.3) 

B  7 

1-6+9 

2.8  (6.4) 

0.95  (2.0) 

2.8  (6.0) 

C  8 

1-7+8 

3.4  (5.0) 

1.46(1.7) 

4.5  t5.6) 

D  8 

1-7+10 

6.8  (6.8) 

1.74  (1.8) 

5.7  (6.0) 

Energy  of  Formation  [mRy/atom] 


Figure  2.14  Calculated  formation  energies  for  selected  compounds 
and  random  bcc  and  hep  alloys  in  Zr-Ru. 
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The  cluster  expansion  coefficients  of  the  excess  volume  and 
bulk  moduli  in  Zr-Ru  are  shown  in  Table  2.8.  The  interactions 
converge  rapidly  beyond  pair  clusters  as  with  Zr-Nb  and  Nb-Ru.  The 
value  of  the  empty  cluster  interaction  for  the  bulk  modulus,  -0.087 
Mbar  is  the  excess  bulk  modulus  of  the  random  alloy  with  an  equi- 
atomic  concentration.  Thus,  for  a  given  lattice  type  it  appears  that 
the  effect  of  order  on  the  bulk  modulus  is  small,  but  the  relaxed  L10 
structure  has  a  bulk  modulus  of  2.633  Mbar,  0.2  Mbar  higher  than 
the  random  bcc  alloy  modulus  --  a  significant  difference. 


Table  2.8.  The  Interaction  coefficients  for  cluster  expansion  of  bulk  modulus  and 
the  volume  of  (a)  bcc-  and  (b)  hep-based  structures  In  Zr-Ru.  The  cluster  labels  are 
from  Figure  2.1. 
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Cluster 

Effective  Interactions 

(label) 

Bj  [mbar] 

Vj  (a.u.3) 

Empty  (1) 

-0.087120 

-5.228230 

Point  (2) 

0.054500 

-1.210090 

Pair  (3) 

-0.002370 

0.797773 

2nd  pair  (4) 

0.026583 

0.743473 

Triangle  (5) 

-0.004540 

0.100842 

Tetrahedron  (6) 

0.002812 

-0.032210 

(b) 


Cluster 

Effective  Interactions 

(label) 

Bi  [mbar] 

Vj  [a.u.3] 

Empty  (1) 

-0.266750 

-3.786960 

Point  (2) 

0.049500 

-0.627460 

Pair  (3  4) 

0.041750 

0.656809 

Triangle  (6  7  8) 

-0.006190 

0.078431 

Tetrahedron  (9) 

0.008125 

-0.076950 
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As  seen  from  this  section  the  effective  interactions  extracted 
from  the  Morse-type  equations  of  state  via  Equation  2.8  provide  an 
adequate  description  of  the  alloy  energy.  Clearly,  a  more 
comprehensive  treatment  of  relaxation  is  required  to  advance  the 
accuracy  of  the  cluster  expansion  to  complete  reliability.  Additionally 
the  use  of  a  full-potential  LDA  total- energy  calculation  is  needed  in 
the  cases  where  structural  relaxation  is  significant;  especially  for 
hexagonal  metals.  The  optimization  criteria  for  selecting  the 
expansion  basis  was  indispensable.  For  even  these  small  sets  of 
structures  most  combinations  of  basis  sets  produced  inadequate 
expansions.  There  remains  no  rigorous  selection  method  for  the  basis 
structures  or  clusters.  Yet.  combining  the  simple  criterion  of  the 
maximum  determinant  of  the  correlation  matrix  in  combination  with 
expanding  the  size  of  the  basis  by  selecting  clusters  in  an.  as  yet.  non¬ 
rigorous  heirarchy  of  increasing  cluster  range  and  complexity  does 
produce  sufficiently  accurate  expansion  for  phase  stability 
calculations. 


CHAPTER  3.  FINITE  TEMPERATURE  CALCULATIONS 


The  solution  to  the  Ising  model  of  the  alloy  free  energy  is 
accomplished  using  the  CVM.  Section  3. 1  outlines  the  assumptions 
used  to  incorporate  vibrational  free  energy  into  the  CVM.  The 
expression  for  configurational  entropy  is  developed  in  Section  3.2. 
The  solution  technique  is  outlined  in  Section  3.3.  and  calculated 
phase  equilibrium  diagram  for  Zr-Nb  is  compared  with  reported 
measurements. 


3,1  Vibrational  free  energy 

In  the  Debye  approximation,  the  velocity  of  sound  is 
constant  for  each  polarization  as  in  the  case  of  a  classical  elastic 
medium.  The  resulting  linear  dispersion  relation  leads  to  the  well- 
known  result  (see,  for  example,  Kittel  [40]): 


cod  =  (6n2N/V)l/3v  , 


[3.1] 


for  a  crystal  with  N  primitive  cells  of  volume  V,  and  cod  is  the 
characteristic  cut-off  frequency.  Further,  if  the  velocity  is  assumed  to 
be  independent  of  the  polarization,  the  characteristic  frequency  is 
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related  to  a  single  characteristic  temperature,  the  Debye 
temperature,  by: 

h  cod  *  k  ©□  ■  (3.2] 

where  h  is  Planck's  constant,  h,  divided  by  2k,  and  k  is  Boltzmann's 
constant.  Using  the  the  velocity  of  sound  in  an  elastic  isotropic 
continuum: 


v=[B/pp/2. 


[3.3] 


where  B  is  the  bulk  modulus  and  p  the  density  leads  to: 


Sfl  If)  134  31 

/  _\I/2 

=  67.48  TO  ,  [3.4  b] 

for  the  Wlgner-Seitz  radius  r  in  atomic  units,  the  atomic  mass,  M  in 
kg.  and  the  bulk  modulus  in  kbar,  the  coefficient  67.48  has  units  of 
Kelvins  -  seconds.  Equation  3.4  b  leads  to  estimates  of  Debye 
temperatures  much  larger  than  experimentally  determined  values. 
The  discrepancy  arises  from  the  treatment  of  the  velocity  of  sound. 
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Anderson  [66]  writes  an  average  velocity  of  sound  as: 


v  = 


.2_  + J_ 

v3  V,3 


m-i/3 


[3.5  a] 


where. 


vt  =  V57p. 


[3.5  b] 


and 


v,  =  Vt7p,  [3.5  c] 

for  longitudinal  and  shear  moduli,  L  and  S.  Moruzzi  et  al  [35] 
combine  Equations  3.5  a-c  with  the  observation  of  Anderson  that  for 
the  non-magnetic  cubic  elements: 

L  *  1.42  B  .  [3.6  a] 


which  is  equivalent  to  the  assumption  that  Poisson's  ratio  is 
approximately  1/3.  Figure  3.1  a  shows  longitudinal  and  bulk  moduli 
for  several  cubic  transition  metals  along  with  the  functional  result  of 
Equation  3.6.  Anderson's  data  also  imply: 


S  =  0.30  B  , 


[3.6  b] 
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(see  Figure  3. 1  b).  From  these  two  assumptions  an  average  velocity  of 
sound  is: 


v  =  0.617 


(3.7] 


which  leads  upon  substitution  into  Equation  3.4  to: 

i  n\,/2 

e“'41-63(M)  • 


[3.8] 


Figure  3.2  compares  the  values  predicted  by  Equation  3.8  with 
experimentally  determined  Debye  temperature  for  the  cubic  non¬ 
magnetic  elemental  metals.  Notably,  a  least-squares  fit  through  the 
measured  values  yields  an  effective  coefficient  of  46.23  K-sec. 

The  anharmonicity  of  the  binding  energy  curve  leads  to 
volume-dependent  effects  in  the  vibration  spectrum  that  can  be 
described  simply  by  the  Gruneisen  relation: 


0  In  Op 

T"  9 In V  * 


or 


[3.10  a] 


a>  =  (a>)o[^f 


[3.10  b] 
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Figure  3.1  a.  Longitudinal  versus  bulk  modulus  for  non¬ 
magnetic  cubic  metals  demonstrating  the  validity  of 
Equation  3.6  a.  (From  Refence  35). 


Figure  3.1  b.  Shear  versus  bulk  modulus  for  non-magnetic 
cubic  metals  compared  with  relation  of  Equation  3.6  b. 
(From  Reference  35). 


Figure  3.2.  Comparison  of  the  measured  values  of  Debye 
temperatures  with  those  predicted  by  Equation  3.8  for  the 
cubic  non-magnetic  transition  metals. 
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If  a  constant  (volume-independent)  Poisson’s  ratio  is  assumed,  then: 


1  i  d  In  B 
6  ’  2  ain  V  * 


[3.11  a] 


or: 


£  V  o2P/9  V2 

3  2g  p/av 


[3.11  b) 


The  validity  of  this  relation  rests  upon  the  assumption  that  all 
vibrational  modes  are  equally  excited  which  is,  effectively,  a  high- 
temperature  average.  Moruzzi  et  al  use  a  low-temperature  average  of 
the  Gruneisen  constant  on  the  grounds  that  most  thermal  expansion 
occurs  at  low  temperatures.  Accordingly,  they  refer  to  the  result  of 
Barron  [67[,  that  by  defining  appropriate  high-  and  low-temperature 
averages  of  the  Gruneisen  constant,  yht  and  y  it: 

Ykt  ■  Ylt  “  g-  •  [3.12] 


Combining  Equation  3.12  with  Equation  3. 1 1  b: 


.  vd2P/dV2 
Y  =  Ylt  =  -  1-  f  . 

2  d  P/9V 


[3.13) 


Ill 


On  the  same  grounds  as  Moruzzi  et  al,  Equation  3.13  is  used  in  this 
study.  The  calculated  Debye  temperature  for  Nbis343K,  and  the 
GrOneisen  parameter  is  1.65  The  experimentally  determined  values 
are  241  K  and  1.74. 

The  volume  dependence  of  the  Debye  temperature  allows 
incorporation  of  thermal  expansion  into  the  vibrational  free  energy 
(see  for  example.  Kittel  [40]): 


where  xq  =  ©d/T. 

Combining  Equation  3.14  with  the  Morse  potential  allows 
computation  of  equilibrium  volumes  at  finite  temperatures  and. 
therefore,  coefficients  of  thermal  expansion.  Figure  3.3  shows  the 
calculated  thermal  expansion  coefficients  of  Nb  as  a  function  of 
temperature  along  with  measured  values.  Also  shown  are  coefficients 
calculated  by  Moruzzi  et  al. 
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Temperature  [K] 


Figure  3.3.  Calculated  coefficients  of  thermal  expansion  for  bcc  Nb. 

The  solid  line  uses  the  total  energy  results  of  this  study,  and  the  dashed 
line  uses  those  of  Reference  35.  The  measured  values  are  from 
Reference  68. 


The  development  of  Equation  3.8  is  justified  only  for  cubic 
metals.  As  yet.  no  comparably  simple  treatment  for  hexagonal  metals 
is  reported.  Figure  3.4  shows  experimentally  determined  Debye 

(  \l/2 

for  the  hexagonal  transition  metals.  A 

least-squares  fit  yields  an  effective  coefficient  of  57.54  K-sec. 

An  ad  hoc  treatment  of  vibrational  modes  similar  to 
Equation  3.8  for  hexagonal  alloys  is  used  in  this  study.  Using 
Equations  3.8  and  3.14  the  transition  of  pure  zirconium  from  hep  to 
bcc  at  866°C  is  reproduced  by  introducing  an  effective  coefficient  of 
45.67  K-sec.  The  induced  transition  temperature  is  quite  sensitive  to 
the  choice  of  the  coefficient  as  is  shown  in  Figure  3.5.  This  coefficient 
is  adopted  in  this  study  for  simplicity,  pending  a  more  fundamental 
analysis.  The  calculated  thermal  expansion  of  hep  zirconium  is 
shown  in  Figure  3.6.  The  calculated  Debye  temperatures  and 
Gruneisen  parameters  for  each  compound  are  shown  in  Table  3. 1  a-c. 
Tsukamoto  {50]  determined  the  Debye  temperature  of  the  partially 
ordered  Llo  compound  Nbs4Ru46.  The  measured  value  is  370  K.  and 
the  value  computed  for  NbRu  Llo  in  the  present  study  is  370.42  K 
(see  Table  3. 1  a). 


Figure  3.4.  Comparison  of  the  measured  values  of  Debye 
temperatures  of  the  hexagonal  transition  metals  with  the  values 
obtained  by  adjusting  the  coefficient  in  Equation  3.8  to  reproduce, 
a  transition  temperature  of  866°C. 


Transition  temperature  (Celsius) 


280  290  300  310  320 


Debye  temperature  (Kelvin)  of  Zr  -  hep 

Figure  3.5.  Zr  hcp-bcc  transition  temperature  as  a  function  of  the 
Debye  temperature  (90)  of  the  hep  phase  with  ©0  of  the  bee  phase 

=  268.93  K.  [The  dotted  line  indicates  the  measured  transition 
temperature,  866°C]. 


Coefficient  of  thermal  expansion  [1  .e-6]  Coefficient  of  thermal  expansion  [1  .e-6] 


Figure  3.6.  Calculated  coefficients  of  thermal  expansion  of  pure  Zr. 
The  top  graph  is  for  the  hep  phase.  The  bottom  graph  shows  the 
behavior  of  both  the  hep  and  bcc  phases. 


Table  3. 1  a.  Calculated  Debye  temperatures  and  Gruneisen  parameters  for  each 
compound  In  Nb-Ru  studied. 
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Compound/ 

Structure 

Debye 

temperature 

[K] 

GrOneisen 

parameter 

Nb/hcp 

369.52 

1.66 

Nb3Ru/DOi9 

385.86 

1.72 

NbRu/Pmma 

405.10 

1.80 

NbRu3/DOig 

424.29 

1.89 

Ru/hcp 

443.94 

2.01 

NbRu/LiRh 

400.64 

1.78 

NbRu/Pmmn 

401.59 

1.79 

NbRu/Hl 

403.77 

1.79 

Nb/fcc 

337.98 

1.67 

Nb3Ru/Ll2 

352.38 

1.73 

NbRu/Llo 

370.42 

1.80 

NbRu3/Ll2 

387.40 

1.89 

Ru/fcc 

402.74 

2.01 

Nb3Ru/D022 

352.19 

1.73 

NbRU3/D022 

387.05 

1.90 

NbRu/NbP 

371.00 

1.81 

Table  3.1a  continued. 


Compound/ 

Structure 

Debye 

temperature 

(K1 

GrOneisen 

parameter 

Nb/bcc 

343.43 

1.65 

Nb3Ru/D03 

360.29 

1.72 

NbRu/B2 

373.92 

1.81 

NbRu/B32 

368.20 

1.79 

NbRu3/D03 

379.62 

1.90 

Ru/bcc 

389.82 

2.02 

NbRu/Tl 

370.11 

1.80 

NbRu/T2 

372.64 

1.82 

NbRu/T3 

372.64 

1.82 

NbRu/T4 

368.20 

1.80 

Nb3Ru/T5 

359.74 

1.73 

NbRu3/T5 

380.43 

1.90 

Table  3. 1  b.  Calculated  Debye  temperatures  and  Gruneisen  parameters  for  each 
compound  in  Zr-Ru  studied. 
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Compound/ 

Structure 

Debye 

temperature 

IK) 

GrOneisen 

parameter 

Zr/hcp 

300.07 

1.46 

Zr3Ru/DOi9 

326.08 

1.55 

ZrRu/Pmma 

359.44 

1.66 

ZrRu3/DOi9 

394.35 

1.80 

Ru/hcp 

443.94 

2.01 

ZrRu/LiRh 

356.76 

1.63 

ZrRu/Pmmn 

352.16 

1.64 

ZrRu/Hl 

354.78 

1.66 

Zr/bcc 

268.95 

1.43 

Zr 3RU/DO3 

305.32 

1.56 

ZrRu/B2 

344.17 

1.70 

ZrRu/B32 

329.49 

1.65 

ZrRu3/D03 

357.18 

1.81 

Ru/bcc 

389.82 

2.02 

ZrRu/Tl 

327.91 

1.65 

ZrRu/T2 

328.16 

1.66 

ZrRu/T3 

325.50 

1.65 

ZrRu/T4 

321.03 

1.63 

Zr3Ru/T5 

299.82 

1.54 

ZrRu3/T5 

359.41 

1.82 

'Table  3.1  c.  Calculated  Debye  temperatures  and  Gruneisen  parameters  for  selected 
structures  In  Zr-Nb. 
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Compound/ 

Structure 

Debye 

temperature 

|K1 

GrOnelsen 

parameter 

Zr/hcp 

301.07 

1.46 

Zr3Nb/DOig 

313.12 

1.49 

ZrNb/Pmma 

329.82 

1.55 

ZrNb3/DOig 

351.60 

1.55 

Nb/hcp 

372.88 

1.60 

ZrNb/LlRh 

325.07 

1.66 

ZrNb/Pmmn 

330.55 

1.56 

Zr/bcc 

268.93 

1.43 

Zr3Nb/D03 

285.43 

1.48 

ZrNb/B2 

302.24 

1.54 

ZrNb/B32 

304.24 

1.54 

ZrNb3/D03 

322.82 

1.60 

Nb/bcc 

343.44 

1.65 

ZrNb/Tl 

303.48 

1.54 

ZrNb/T2 

302.61 

1.54 

ZrNb/T3 

302.77 

1.53 
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In  this  study  the  configurational  entropy  is  calculated  using 
the  cluster  variation  method  (CVM).  The  CVM  is  based  on  the  fact 
that  for  a  given  probability  distribution  X(o),  the  exact  configurational 
entropy  is  given  by  (69): 

S=  -kB  X  X(a)lnX(a) .  (3.15) 

o 

where  the  sum  is  over  all  2N  configurations  of  the  crystal.  A  sequence 
of  cluster  entropies  defined  by: 

S«=  -kB  X  Xa(aa)lnXa(aa) ,  [3.16] 

a 

clearly  converges  to  the  exact  configurational  entropy  as  the  cluster 
size  a  approaches  N.  Sanchez  and  de  Fontaine  [69]  used  a  MObius 
transformation  to  write  the  cluster  entropies,  Soc.  in  terms  of  a  set  of 
irreducible  cluster  contributions,  Sp: 

So=  X  Sp  * 

Pea 


[3.17] 


where  the  sum  extends  over  all  subclusters  of  a,  including  a  but 
excluding  the  empty  cluster.  As  in  all  cluster  expansions  the  key 
approximation  of  the  CVM  is  to  neglect  the  contributions  Sa  for 
clusters  larger  than  a  maximum  cluster.  Due  to  the  symmetry  of  the 
space  group  of  the  crystal,  Equation  3. 15  can  be  written  as: 

m 

S=£z„Sn.  [3.18] 

n=l 

where  m  labels  the  maximum  cluster,  and  zn  is  the  number  of  n-type 
clusters  per  lattice  site.  This  equation,  in  turn,  may  be  rewritten  in 
terms  of  cluster  entropies: 

m  m 

S  -N  £  zn  an  Sn  =  -Nka  £  zn  an£xn(on)  In  Xn(cn).  [3.19] 

n=l  n=l  on 

where  the  coefficients  an  are  determined  using: 

£  ap  =  1  ,  [3.20] 

Pea 

which  is  valid  for  each  subcluster  a  of  the  maximum  cluster.  In 
Equation  3.20  the  sum  extends  over  all  subclusters  p  of  the  maximum 
clusters  that  contain  a. 
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3.3  CVM  solution  method 

Combining  Equations  2.5.  3.14.  and  3.19  into  an  expression 
for  the  alloy  free  energy  functional  [14): 

=  f  Zn  JnfV.T) 

n-0 

+  kT  I  Z„  an  I  XJoJ  In  XJoJ  .  [3.211 

n»  1  On 


where  the  cluster  interactions  Jn  are  obtained  from  a  cluster 
expansion  of  the  excess  structural  energy  and  the  excess  vibrational 
free  energy.  The  alloy  free  energy  expression  depends  explicitly  only 
upon  the  volume  V.  the  temperature,  and  the  correlations  £n.  where 
the  cluster  probabilities  Xn  may  be  determined  by  Equation  2.9: 


Xn(o)  =  ± 
2 


1  +  X  Anp(o) 
P-I 


[3.22] 


Then,  for  a  given  temperature  and  concentration,  the  equilibrium  free 
energy  is  obtained  by  minimizing  the  free  energy  functional  with 
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respect  to  the  volume  and  the  correlations.  Recall  that  the  point 
correlation. 

4i  =  2xj  -1  =  xj-  X2  (2.22) 

is  fixed  for  a  given  concentration  xi  (=  c). 

Calculation  of  phase  boundaries  is  accomplished  efficiently 
using  the  grand  potential  scheme  proposed  by  Kikuchi  [70] .  The 
grand  potential  is  obtained  by  the  Legendre  transformation  of  the  free 
energy: 


Q(H.  T.  U),  V)  =  F(T,  {  4},V)+  (i£j  ,  [3.23] 

where  {  £}  is  the  set  of  all  correlations,  and  n  is  the  effective  chemical 
potential: 


V-  =  - 


9F 


[3.24] 


since  by  definition: 

F  =  p,xi  +  p2x2 


[3.25] 


and. 
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3F  _  9F^  3xj_  9F  dx2 

35, '  a*i  35,  35, 

=  (Hi  -  M-i)  M-  13.26] 


Similarly, 


Q  =  F  -  |i(x,  -  xj 


Cixi  +  N 


[3.27] 


Therefore  the  conditions  for  equilibrium  between  phases  a  and  p  is 
that: 


Q«  =  QP 


[3.28] 


and 


pa  =  |iP  .  [3.29] 

The  solution  procedure  used  most  often  for  determining 
phase  boundaries  in  this  study  is  to  fix  temperature  and  effective 
chemical  potential  fi,  then  using  a  Newton-  Raphson  iterative 
algorithm  to  determine  the  equilibrium  variables  V  and  {  ^}.  This 
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algorithm  Includes  an  inner  loop  In  which  Qa  and  QP  each  are 
minimized  using  a  joint  minimization  scheme.  That  is,  for  each  phase 
the  minimization: 


an 

aiv.  $ 


=o  . 


[3.30] 


is  achieved  using  the  Newton-Raphson  correction: 


[V.  ynew  =  {V,  $)°'d  -  H'1 


an 


an-\  $} 


[3.31] 


where  the  Hessian  H  is: 

a2n  a2n 

H=  as, 

a2n  a2n 
a^,  av  av2 


[3.14] 


The  outer  loop  iterates  on  (i  until  Equations  3.28  and  3.29  are 
satisfied  self-consistently. 

The  bcc  miscibility  gap  in  Zr-Nb  was  calculated  using  an 
Ising  Hamiltonian  varying  in  maximum  interaction  range  from  second 
nearest  neighbors  to  fifth  nearest  neighbors  and  the  CVM 
approximation  with  maximum  cluster  size  ranging  from  the 
tetrahedron  to  the  joint  9  point-6  point  cluster.  The  results,  shown  in 
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Figure  3.7  indicate  rapid  convergence  with  the  interaction  r singe. 
Curves  1  through  4  are  calculated  using  global  relaxation.  Curve  1 
denotes  the  lowest  order  approximation  of  both  the  Hamiltonian  and 
the  statistical  range  —  2nd  nearest  neighbors,  i.  e.,  cluster  set  A  as 
defined  in  Table  2.2.  If  the  statistical  range  is  increased  to  include  the 
Joint  9  point  -  6  point  cluster,  but  the  Hamiltonian  is  confined  to  the 
same  basis  set  as  Curve  1,  there  is  very  little  change  --  as  seen  by 
comparing  curves  1  and  2.  If  the  Hamiltonian  range  is  increased  to 
include  third  nearest  neighbors  (basis  set  B  in  Table  2.2),  there  is  a 
large  decrease  in  the  critical  temperature  of  the  miscibility  gap,  curve 
4.  Inclusion  of  fifth  nearest-neighbor  interactions  (basis  set  C  Table 
2.2)  produces  little  change  (curve  3)  in  the  calculated  miscibility  gap. 
Thus,  the  CVM  calculation  of  the  miscibility  gap  is  deemed  converged 
at  3rd  nearest  neighbors.  Using  the  same  interaction  range  the 
miscibility  gap  was  calculated  using  totally  relaxed  interactions.  The 
result,  curve  5.  is  significantly  different  from  the  well  converged 
globally  relaxed  curve  4.  A  more  physically  representative  relaxation 
method  should  lie  somewhere  between  the  total  and  global  relaxation 
methods,  but  it  should  be  reiterated  that  the  globally  relaxed  cluster 
expansion  reproduced  total  energies  very  precisely. 

The  equilibrium  of  the  hexagonal  phase  was  calculated  using 
the  Joint  tetrahedron-  octahedron  cluster  with  second  nearest 
neighbor  interactions  included.  As  discussed  in  Section  3.1  the 
Debye  temperatures  of  hexagonal  alloys  have  been  determined  using  a 
coefficient  in  Equation  3.8  which  reproduces  the  measured  transition 
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temperature  of  pure  Zr.  This  treatment  yielded  the  phase  diagram 
shown  in  Figure  3.8.  Also  shown  are  several  experimental  data. 
Global  relaxation  is  assumed  in  the  calculation.  The  computed 
critical  temperature  of  the  miscibility  gap  is  975  °C  at  62%  Nb 
compared  with  the  reported  values  of  988  °C  and  61%  Nb  [7].  The 
computed  monotectoid  equilibrium  temperature  is  855  °C  compared 
with  the  experimental  893  °C. 


Figure  3.7.  Zr-Nb  miscibility  gap  calculated  using  different 
approximations  of  the  CVM: 

1  -  2nd  nearest  neighbors;  tetrahedron 

2  -  2nd  nearest  neighbors;  9  point-  6  point 

3  -  5*  nearest  neighbors;  9  point  -  6  point 

4  -  3rd  nearest  neighbors;  9  point  -  6  point 

5  -  3rd  nearest  neighbors;  9  point  -  6  point;  totally  relaxed 


No 


Temperature  (K) 
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Concentration  of  Niobium 


Figure  3.8.  Calculated  phase  diagram  for  Zr-Nb  along  with 
data  from  References  41  and  42.  The  calculations  include  third 
nearest  neighbors  in  the  cluster  expansion  and  globally 
relaxed  interactions. 


CONCLUSIONS 


The  phase  stability  of  alloys  in  each  of  the  three  systems  Zr-Nb, 
Nb-Ru,  and  Ru-Zr  was  examined  using  the  cluster  expansion 
technique  described  in  Chapter  2.  The  equation  of  state  (Equation 
1.5)  for  each  of  several  perfectly  ordered  structures  in  each  system 
was  computed  using  the  local-density  approximation  (2)  to  density- 
functional  theoiy  (I).  The  relatively  short-range  effective  chemical 
interactions  were  determined  using  the  inversion  of  the  cluster 
expansion  (Equation  2.9)  of  a  set  of  total  energies  which,  in  general, 
may  depend  upon  very  long-range  interactions.  The  convergence  of 
the  expansion  and  the  sensitivity  of  the  inversion  to  errors  are 
coupled  problems  of  selecting  optimum  basis  sets  of  clusters  and 
structures  and  are  discussed  in  Sections  2.3-2.6.  An  efficient  scheme 
for  obtaining  a  sufficiently  converged  expansion  was  developed  and 
applied  to  the  three  alloy  systems.  The  method  consists  of  an 
arbitrary  initial  choice  for  the  maximum  interaction  range  for  a  given 
Bravais  lattice.  For  this  maximum  range  there  is  a  unique  maximum 
cluster  (or  set  of  joint  clusters)  from  which  a  set  of  sub-clusters  are 
chosen  to  form  the  cluster  basis.  An  optimum  choice  of  the  set  of 
basis  structures  is  then  determined  using  the  maximum  determinant 
method,  and  the  total  energies  are  calculated  using  the  local  density 
approximation.  For  example,  in  the  case  of  structures  with  the  hep 
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parent  lattice  a  maximum  range  of  nearest  neighbors  was  first 
chosen,  for  which,  the  associated  largest  cluster  is  a  regular 
tetrahedron.  The  eight  different  sub-clusters  of  this  tetrahedron  were 
reduced  to  five  by  assuming  the  isotropy  (or  equivalence)  of  the 
interactions  among  all  three  types  of  triangles  and  between  the  two 
types  of  pairs.  Since  the  vertices  of  the  configurational  polyhedron  for 
this  choice  of  cluster  set  is  known,  the  choice  of  basis  structures  was 
obvious:  the  structures  corresponding  to  those  vertices.  However,  for 
Zr-Nb  the  interactions  derived  using  this  representation  were 
insufficient  to  reproduce  the  formation  energy  of  the  test  structure  so 
the  cluster  basis  was  expanded  to  include  second-nearest  neighbors 
with  the  regular  tetrahedron  retained  as  the  largest  cluster  size.  The 
corresponding  structure  set  was  determined  using  the  maximum 
determinant  method  (Section  2.5),  and  the  resulting  cluster  expansion 
was  deemed  adequate  by  the  convergence  criterion  (Section  2,6). 

In  this  study,  the  initial  choices  for  and  subsequent 
augmentations  of  the  basis  cluster  sets  were  made  by  somewhat  ad 
hoc  methods.  A  fully  rigorous  implementation  would  consist  of 
exhaustive  searches  of  each  combination  of  sub-clusters  for  a  given 
interaction  range  and  the  use  of  large-scale  vertex  enumeration 
studies  neither  of  which  are  within  the  scope  of  this  study.  The  three 
criteria  for  choosing  the  basis  structures  yielded  nearly  identically 
converged  expansions.  As  the  determinant  method  is  the  simplest  of 
the  three,  it  was  used  as  the  criterion  in  this  study.  The  results 
presented  in  Section  2.  6  indicate  the  precision  of  the  cluster 
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expansions  for  these  alloy  systems  as  better  than  5%  of  the  formation 
energy. 

Using  the  converged  cluster  expansions,  limited  ground-state 
searches  were  performed  by  calculating  the  formation  energy  of  every 
possible  ordered  structure  for  an  a  priori  choice  of  maximum  unit  cell 
size.  Using  this  method,  the  ground  states  of  each  of  the  three 
systems  were  determined.  Notably,  the  ground-state  structure  of 
NbRu3  was  predicted  to  be  the  DO  19  structure.  The  calculated  total 
energies  of  the  DO19  and  LI2  structures  are  very  close,  and  since  no 
relaxation  of  the  c/a  parameter  of  the  hexagonal  structure  is  allowed 
within  the  atomic-sphere  approximation,  a  full-potential  total  energy 
calculation  is  necessary  for  calculation  of  the  margin  of  stability 
between  the  cubic  and  hexagonal  structures. 

An  additional  stable  ground-state  phase.  Nb3Ru  with  the  DO3 
crystal  structure,  is  predicted  using  the  cluster  expansion.  There  is 
no  experimental  evidence  of  an  ordered  phase  at  this  composition,  but 
this  may  due  to  the  very  similar  scattering  factors  of  Nb  and  Ru  — 
which  makes  experimental  determination  of  the  ordering  of  the  other 
phases.  NbRu  (B2  and  Llo)  and  NbRu3  difficult  as  well. 

The  cluster  variation  method  is  outlined  in  Sections  3.2  and  3.3 
for  calculation  of  the  configurational  free  energy  at  finite 
temperatures.  Figure  3.7  indicates  the  rapid  convergence  of  the 
configurational  entropy  with  the  tetrahedron-octahedron  or  the  nine 
point-six  point  approximations  for  the  hexagonal  and  bcc  phases 
respectively. 
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The  contribution  of  the  vibrational  entropy  to  the  excess  free 
energy  is  usually  insignificant  for  alloys  of  the  same  lattice  type. 
However,  for  phase  transformations  in  which  the  lattice  type  changes, 
the  vibrational  free  energy  may  be  significant.  This  is  the  case  for 
pure  Zr,  which  undergoes  a  transition  from  hep  to  bcc  at  866°C.  In 
this  study,  vibrational  free  energy  is  calculated  using  the  Debye- 
Gruneisen  approximation,  and  Debye  temperatures  calculated  from 
an  effective  speed  of  sound  following  Moruzzi  et  al  [35].  This 
treatment  allows  the  calculation  of  coefficients  of  thermal  expansion 
which  for  niobium  compare  very  well  with  measurements  (Figure  3.3). 
The  utility  of  this  approach  rests  upon  incorporating  an  effective 
speed  of  sound  for  the  phase.  Moruzzi  et  al  obtained  the  relation 
(Equation  3.7)  for  cubic  metals  from  data  from  Anderson  [66].  For 
hexagonal  phases  this  treatment  is  not  applicable  so,  in  this  study,  an 
ad  hoc  treatment  was  adopted  which  requires  a  single  adjustable 
parameter.  For  simplicity,  the  parameter  was  adjusted  to  reproduce 
the  actual  transition  temperature  of  Zr. 

The  incorporation  of  vibrational  free  energy  allows  the 
calculation  of  the  full  Zr-Nb  phase  diagram  (Figure  3.8).  The  Ising 
Hamiltonian  for  this  system  is  well  converged  for  both  hep  and  bcc 
alloys  as  is  the  statistical  approximation.  The  computed  equilibria  are 
in  excellent  agreement  with  reported  data:  the  critical  point  of  the 
miscibility  gap  is  975°C  at  62%  niobium,  the  measured  values  988°C 
and  61%,  respectively. 
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